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: user-unsubscr...@commons.apache.org
For additional commands, e-mail: user-h...@commons.apache.org
