yanyankangkang / efficient-java-matrix-library

Automatically exported from code.google.com/p/efficient-java-matrix-library
0 stars 0 forks source link

Accuracy loss #24

Closed GoogleCodeExporter closed 9 years ago

GoogleCodeExporter commented 9 years ago
Hi!
I'am using EJML for my numerical optimization software. I think there is some 
problem with accuracy in calculations (i.e. innerProduct)

What steps will reproduce the problem?

public class PrecisionTest
{
  public static void main(String[] arg0)
  {
    int testCount= 10;

    for(int testIndex=0; testIndex<testCount; testIndex++) 
    {
      double s1= 0.0, s2= 0.0;

      double[] p1= new double[50], p2= new double[50];

      Arrays.fill(p1,Math.PI); Arrays.fill(p2,Math.PI);

      for(int testRepeat=0; testRepeat<10000; testRepeat++)
      {
    for(int i=0; i<50; i++)
    {
      s2 += p1[i]*p2[i];
    }
      }

      DenseMatrix64F v1= new DenseMatrix64F(1,50,true,p1), v2= new DenseMatrix64F(50,1,false,p2);

      for(int testRepeat=0; testRepeat<10000; testRepeat++)
      {
    s1 += VectorVectorMult.innerProd(v1, v2);
      }

      System.out.println("s1="+String.format("%.25f",s1)
                         +",s2="+String.format("%.25f",s2)
                         +",s1-s2="+String.format("%.25f",s1-s2));
    }
  }
}

What is the expected output? What do you see instead?

Expected Abs(s1-s2)= 0.0, but realy Abs(s1-s2) > 0.0. ??? please check the loss 
of accuracy in calculations. I'll be happy to be wrong in its conclusions.

What version of the product are you using? On what operating system?

EJML-0.18.jar, source: maven repo.

Please provide any additional information below.

Original issue reported on code.google.com by Azimu...@gmail.com on 12 Feb 2012 at 8:33

GoogleCodeExporter commented 9 years ago
Windows 7 64-bit, CPU: AMD Athlon X2

Original comment by Azimu...@gmail.com on 12 Feb 2012 at 8:46

GoogleCodeExporter commented 9 years ago
No bug, finite precision math is tricky and highly dependent on the order in 
which operations happen.  The two loops are "mathematically" equivalent but not 
computationally equivalent.

To get the same solution in both loops simply sum up the batch of 50 first, 
then add it to s2:

double sum = 0;
for(int i=0; i<50; i++)
{
    sum += p1[i]*p2[i];
}s2 += sum;

Original comment by peter.ab...@gmail.com on 14 Feb 2012 at 8:56

GoogleCodeExporter commented 9 years ago
thanks

Original comment by Azimu...@gmail.com on 15 Feb 2012 at 10:37