Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Double precision inner product

Subject: Double precision inner product

From: Matt J

Date: 17 Jan, 2011 16:29:07

Message: 1 of 22

I have two large arrays of type single, X and Y, and I need their inner product
X(:).'*Y(:) calculated to double precision.

I cannot see a way of doing this without first casting X and Y to double, which because of their large size, I do not want to do. The closest that MATLAB seems to offer is

sum(X(:).*Y(:),'double')

which is not ideal since
(a) the multiplies are still done in single precision
(b) It creates a large intermediate array X(:).*Y(:)

Does anyone know of a FEX tool that will accomplish this?

Subject: Double precision inner product

From: Bruno Luong

Date: 17 Jan, 2011 17:57:05

Message: 2 of 22

"Matt J" wrote in message <ih1qoi$ko6$1@fred.mathworks.com>...

>
> Does anyone know of a FEX tool that will accomplish this?

No need FEX tool, it is relatively straightforward to code in mex. Fell free to OPENMP for multithreading sum.

/*************************************************************************
 * MATLAB MEX ROUTINE accusingledot.c
 *
 * > accusingledot(A,B)
 * performs dot product of two single array
 *
 * User must make sure A is not shared by other array.
 * B must have the same size as A, no check will be carried out
 *
 * A, B must be of class single
 *
 * > mex -largeArrayDims accusingledot.c % on 64 bit platform
 * > mex accusingledot.c % on 64 bit platform
 ************************************************************************/

#include "mex.h"

#define A prhs[0]
#define B prhs[1]

/* Gateway of inplaceprod */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[]) {
    
    size_t i, n;
    double s;
    mxClassID ClassA, ClassB;
    float *a, *b;
    
    if (nrhs != 2)
        mexErrMsgTxt("accusingledot: two arguments required.");
    /* Get the classes of A and B */
    ClassA = mxGetClassID(A);
    ClassB = mxGetClassID(B);
    
    /* n = numel(A) */
    n = mxGetNumberOfElements(A);
    if (n != mxGetNumberOfElements(B))
        mexErrMsgTxt("accusingledot: A & B must have the same number of elements.");
    if (ClassA == mxSINGLE_CLASS || ClassB == mxSINGLE_CLASS) {
        s = 0;
        a = (float*)mxGetData(A);
        b = (float*)mxGetData(B);
        for (i=0; i<n; i++) s += a[i]*b[i];
    } else
        mexErrMsgTxt("accusingledot: A & B must be single type.");
    
    plhs[0] = mxCreateDoubleScalar(s);
    
    return;
} /* of accusingledot */

% Bruno

Subject: Double precision inner product

From: Bruno Luong

Date: 17 Jan, 2011 18:05:07

Message: 3 of 22

Note, you might change:

> for (i=0; i<n; i++) s += a[i]*b[i];

to

> for (i=0; i<n; i++) s += (double)a[i] * (double)b[i];

If you want to cast even earlier (more accurate product).

Bruno

Subject: Double precision inner product

From: Matt J

Date: 17 Jan, 2011 18:18:06

Message: 4 of 22

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih1vth$1f2$1@fred.mathworks.com>...
>
> * > mex -largeArrayDims accusingledot.c % on 64 bit platform
> * > mex accusingledot.c % on 64 bit platform

Thanks, Bruno. Above, you have 2 different compile instructions both for 64-bit. Was one of them supposed to be for 32-bit?

Subject: Double precision inner product

From: Bruno Luong

Date: 17 Jan, 2011 18:38:05

Message: 5 of 22

"Matt J" wrote in message <ih214t$oot$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih1vth$1f2$1@fred.mathworks.com>...
> >
> > * > mex -largeArrayDims accusingledot.c % on 64 bit platform
> > * > mex accusingledot.c % on 64 bit platform
>
> Thanks, Bruno. Above, you have 2 different compile instructions both for 64-bit. Was one of them supposed to be for 32-bit?

Yes. Sorry the command for 32-bit is WITHOUT -largeArrayDims option.

Bruno

Subject: Double precision inner product

From: Derek O'Connor

Date: 19 Jan, 2011 03:54:05

Message: 6 of 22

"Matt J" wrote in message <ih1qoi$ko6$1@fred.mathworks.com>...
> I have two large arrays of type single, X and Y, and I need their inner product
> X(:).'*Y(:) calculated to double precision.
>
> I cannot see a way of doing this without first casting X and Y to double, which because of their large size, I do not want to do. The closest that MATLAB seems to offer is
>
> sum(X(:).*Y(:),'double')
>
> which is not ideal since
> (a) the multiplies are still done in single precision
> (b) It creates a large intermediate array X(:).*Y(:)
>
> Does anyone know of a FEX tool that will accomplish this?
-------------------------------------

I tested 4 methods with n = 10^6

x = rand(n,1,'single');
y = rand(n,1,'single');

rdot = Dot2(x,y); % Rump's Quad Prec Dot

mdot = dot(x,y); % Matlab

ssumxy = 0.0;
for i = 1:n
    ssumxy = ssumxy + x(i)*y(i);
end;

ddsumxy = 0.0;
for i = 1:n
    ddsumxy = ddsumxy + double(x(i))*double(y(i));
end;

The errors relative to Rump's rdot were

mdot : 3.1320999e-007

ssumxy : 5.7160823e-006

ddsumxy : 0.0


What shocked me was that ssumxy was 1500 times slower than mdot.
ddsumxy was 2.6 times slower than mdot.

When will we get REAL single precision?

Derek O'Connor


System: Dell Precision 690, Intel 2xQuad-Core E5345 @ 2.33GHz
        16GB RAM, Windows7 64-bit Professional.
MATLAB Version 7.6.0.324 (R2008a)

Subject: Double precision inner product

From: Bruno Luong

Date: 19 Jan, 2011 08:12:05

Message: 7 of 22

"Derek O'Connor" wrote in message <ih5n8t$ark$1@fred.mathworks.com>...
> "Matt J" wrote in message <ih1qoi$ko6$1@fred.mathworks.com>...
> > I have two large arrays of type single, X and Y, and I need their inner product
> > X(:).'*Y(:) calculated to double precision.
> >
> > I cannot see a way of doing this without first casting X and Y to double, which because of their large size, I do not want to do. The closest that MATLAB seems to offer is
> >
> > sum(X(:).*Y(:),'double')
> >
> > which is not ideal since
> > (a) the multiplies are still done in single precision
> > (b) It creates a large intermediate array X(:).*Y(:)
> >
> > Does anyone know of a FEX tool that will accomplish this?
> -------------------------------------
>
> I tested 4 methods with n = 10^6
>
> x = rand(n,1,'single');
> y = rand(n,1,'single');
>
> rdot = Dot2(x,y); % Rump's Quad Prec Dot
>
> mdot = dot(x,y); % Matlab
>
> ssumxy = 0.0;
> for i = 1:n
> ssumxy = ssumxy + x(i)*y(i);
> end;
>
> ddsumxy = 0.0;
> for i = 1:n
> ddsumxy = ddsumxy + double(x(i))*double(y(i));
> end;
>
> The errors relative to Rump's rdot were
>
> mdot : 3.1320999e-007
>
> ssumxy : 5.7160823e-006
>
> ddsumxy : 0.0
>
>
> What shocked me was that ssumxy was 1500 times slower than mdot.
> ddsumxy was 2.6 times slower than mdot.
>
> When will we get REAL single precision?

The result is different to the test I carried out on 2011A Prerelease:

ddsumxy gives different result than Matlab dot (I limit 'n' to 80000 to avoid pollution by multithreading effect on sum).

Bruno

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 08:30:21

Message: 8 of 22

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih66cl$t5g$1@fred.mathworks.com>...
>
> The result is different to the test I carried out on 2011A Prerelease:
>
> ddsumxy gives different result than Matlab dot (I limit 'n' to 80000 to avoid pollution by multithreading effect on sum).
>
> Bruno

FYI, the MATLAB built-in dot function is rather slow and less accurate than just doing a matrix multiply. i.e., it is better to just to x'*y rather than dot(x,y). Also, the MTIMESX SPEED mode and SPEEDOMP modes can often meet or beat MATLAB for speed and accuracy for a large dot product calculation (MTIMESX has hand-coded unrolled loops).

For the specific case of single inputs but calculation and accumulation in double, the BLAS has a specific routine for doing this calculation, DSDOT, but I cannot find it at the moment in the supplied BLAS library. I will look some more tomorrow.

James Tursa

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 08:50:04

Message: 9 of 22

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih1vth$1f2$1@fred.mathworks.com>...

(snip)
> if (ClassA == mxSINGLE_CLASS || ClassB == mxSINGLE_CLASS) {

The above should be:

     if (ClassA == mxSINGLE_CLASS && ClassB == mxSINGLE_CLASS) {

James Tursa

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 08:54:05

Message: 10 of 22

"James Tursa" wrote in message <ih67et$808$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih66cl$t5g$1@fred.mathworks.com>...
> >
> > The result is different to the test I carried out on 2011A Prerelease:
> >
> > ddsumxy gives different result than Matlab dot (I limit 'n' to 80000 to avoid pollution by multithreading effect on sum).
> >
> > Bruno
>
> FYI, the MATLAB built-in dot function is rather slow and less accurate than just doing a matrix multiply. i.e., it is better to just to x'*y rather than dot(x,y). Also, the MTIMESX SPEED mode and SPEEDOMP modes can often meet or beat MATLAB for speed and accuracy for a large dot product calculation (MTIMESX has hand-coded unrolled loops).

P.S. You may have to use LOOPSOMP mode to force the OpenMP loops code to be used.

James Tursa

Subject: Double precision inner product

From: Derek O'Connor

Date: 19 Jan, 2011 15:58:05

Message: 11 of 22

"James Tursa" wrote in message <ih68rd$b06$1@fred.mathworks.com>...
> "James Tursa" wrote in message <ih67et$808$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih66cl$t5g$1@fred.mathworks.com>...
> > >
> > > The result is different to the test I carried out on 2011A Prerelease:
> > >
> > > ddsumxy gives different result than Matlab dot (I limit 'n' to 80000 to avoid pollution by multithreading effect on sum).
> > >
> > > Bruno
> >
> > FYI, the MATLAB built-in dot function is rather slow and less accurate than just doing a matrix multiply. i.e., it is better to just to x'*y rather than dot(x,y). Also, the MTIMESX SPEED mode and SPEEDOMP modes can often meet or beat MATLAB for speed and accuracy for a large dot product calculation (MTIMESX has hand-coded unrolled loops).
>
> P.S. You may have to use LOOPSOMP mode to force the OpenMP loops code to be used.
>
> James Tursa

-----------------------------------

James,

I have re-run these tests with 1 CPU

R2008b 1 CPU
                         x.y RelErr(Dot2) Time (secs) Time rel. x'*y

Dot2(x,y) 2.4970e+005 0 1.2050e+002 3.6358e+004

dot(x,y) 2.4966e+005 1.5894e-004 1.3064e-002 3.9417e+000

x'*y 2.4970e+005 3.3791e-006 3.3142e-003 1.0000e+000

ssumxy 2.4966e+005 1.5894e-004 9.1586e+000 2.7635e+003

dssumxy 5.0002e+005 1.0025e+000 2.2298e-002 6.7281e+000

ddsumxy 2.4970e+005 0 2.0944e-002 6.3194e+000


This confirms what you say about Matlab's dot(x,y) -- slow and
innacurate (rel err = 1.6e-4)

x'*y is fast but inaccurate (rel err = 3.4e-6)

ssumxy is inaccurate and incredibly slow ~ 3000 times slower!

dssumxy = dssumxy + double(x(i)*y(i)) is awful (rel err > 1)

ddsumxy = ddssumxy + double(x(i))*double(y(i)) is slow
but accurate.

My confidence in Matlab has gone down another notch.

Derek O'Connor.

Subject: Double precision inner product

From: Matt J

Date: 19 Jan, 2011 16:54:04

Message: 12 of 22

"James Tursa" wrote in message <ih68js$pfc$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih1vth$1f2$1@fred.mathworks.com>...
>
> (snip)
> > if (ClassA == mxSINGLE_CLASS || ClassB == mxSINGLE_CLASS) {
>
> The above should be:
>
> if (ClassA == mxSINGLE_CLASS && ClassB == mxSINGLE_CLASS) {
>

Could I avoid this type-checking if I make a and b void pointers and cast appropriately?

       a = (void*)mxGetData(A);
       b = (void*)mxGetData(B);

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 18:32:05

Message: 13 of 22

"Matt J" wrote in message <ih74vc$khl$1@fred.mathworks.com>...
> "James Tursa" wrote in message <ih68js$pfc$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ih1vth$1f2$1@fred.mathworks.com>...
> >
> > (snip)
> > > if (ClassA == mxSINGLE_CLASS || ClassB == mxSINGLE_CLASS) {
> >
> > The above should be:
> >
> > if (ClassA == mxSINGLE_CLASS && ClassB == mxSINGLE_CLASS) {
> >
>
> Could I avoid this type-checking if I make a and b void pointers and cast appropriately?
>
> a = (void*)mxGetData(A);
> b = (void*)mxGetData(B);

No. In the first place, mxGetData already returns a (void *) type so your cast above does nothing. But that aside, changing the pointer type does not change the underlying data ... i.e. when you dereference the pointer there is no mechanism in C to automatically do a type conversion based on the actual underlying data type. To get the conversion to happen you need to explicitly code it, like Bruno did in an earlier post.

James Tursa

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 18:36:05

Message: 14 of 22

"James Tursa" wrote in message <ih67et$808$1@fred.mathworks.com>...
>
> For the specific case of single inputs but calculation and accumulation in double, the BLAS has a specific routine for doing this calculation, DSDOT, but I cannot find it at the moment in the supplied BLAS library. I will look some more tomorrow.

I checked into this and it appears that the DSDOT routine was not added to the BLAS library that MATLAB uses until version R2008a. So one could write a simple mex routine to call this, but it would only work in R2008a or later unless you used a different BLAS library.

James Tursa

Subject: Double precision inner product

From: Bruno Luong

Date: 19 Jan, 2011 19:46:07

Message: 15 of 22

"James Tursa" wrote in message <ih68js$pfc$1@fred.mathworks.com>...
>
>
> The above should be:
>
> if (ClassA == mxSINGLE_CLASS && ClassB == mxSINGLE_CLASS) {
>

True. Thanks for the correction.

Bruno

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 21:45:08

Message: 16 of 22

"Derek O'Connor" wrote in message <ih71md$2e0$1@fred.mathworks.com>...
>
> I have re-run these tests with 1 CPU
>
> R2008b 1 CPU
> x.y RelErr(Dot2) Time (secs) Time rel. x'*y
>
> Dot2(x,y) 2.4970e+005 0 1.2050e+002 3.6358e+004
>
> dot(x,y) 2.4966e+005 1.5894e-004 1.3064e-002 3.9417e+000
>
> x'*y 2.4970e+005 3.3791e-006 3.3142e-003 1.0000e+000
>
> ssumxy 2.4966e+005 1.5894e-004 9.1586e+000 2.7635e+003
>
> dssumxy 5.0002e+005 1.0025e+000 2.2298e-002 6.7281e+000
>
> ddsumxy 2.4970e+005 0 2.0944e-002 6.3194e+000
>
>
> This confirms what you say about Matlab's dot(x,y) -- slow and
> innacurate (rel err = 1.6e-4)
>
> x'*y is fast but inaccurate (rel err = 3.4e-6)
>
> ssumxy is inaccurate and incredibly slow ~ 3000 times slower!
>
> dssumxy = dssumxy + double(x(i)*y(i)) is awful (rel err > 1)
>
> ddsumxy = ddssumxy + double(x(i))*double(y(i)) is slow
> but accurate.
>
> My confidence in Matlab has gone down another notch.
>
> Derek O'Connor.

I think you are comparing apples to oranges for some of these. In particular, labeling the x'*y calculation as inaccurate and the ddsumxy calculation as accurate. The x'y calculation may be doing intermediate calculations/storage in single precision whereas the ddsumxy method is doing all intermediate calculations/storage in double precision. In fact, the x'*y result is pretty close to being accurate to the full single precision mantissa length which is quite acceptable performance. So I don't see any accuracy issue with the x'*y method.

I did extensive testing of dot product accuracy calculations when I coded my MTIMESX FEX submission, comparing the results to 100 decimal digit accurate calculations using the same inputs. I always saw the BLAS dot product routine to be among the best for accuracy, although as I mentioned previously you can sometimes get a slightly more accurate result with custom unrolled loops. The amount of unrolling (10) I eventually decided on was a compromise between speed and accuracy based on a variety of different sized inputs. I never saw the BLAS dot product result to be so bad that I would call it inaccurate.

James Tursa

Subject: Double precision inner product

From: Jan Simon

Date: 19 Jan, 2011 22:08:04

Message: 17 of 22

Dear James,

> I always saw the BLAS dot product routine to be among the best for accuracy, although as I mentioned previously you can sometimes get a slightly more accurate result with custom unrolled loops. The amount of unrolling (10) I eventually decided on was a compromise between speed and accuracy based on a variety of different sized inputs.

I do not think, that the effect of loop unrolling to the accuracy can be determined in general without considering the processed data.

> I never saw the BLAS dot product result to be so bad that I would call it inaccurate.

A short example, Matlab 2009a, Windows:
  x = ones(1, 3);
  y = [1e9, 1e-9, -1e9];
  x * y' % Calls BLAS
>> 1e-9
  sum(x .* y)
>> 0

Whenever the BLAS dot product seems to be not accurate enough, try the XBLAS functions, which are stable and well tested also.

Kind regards, Jan

Subject: Double precision inner product

From: James Tursa

Date: 19 Jan, 2011 23:47:04

Message: 18 of 22

"Jan Simon" wrote in message <ih7nc4$4ma$1@fred.mathworks.com>...
> Dear James,
>
> > I always saw the BLAS dot product routine to be among the best for accuracy, although as I mentioned previously you can sometimes get a slightly more accurate result with custom unrolled loops. The amount of unrolling (10) I eventually decided on was a compromise between speed and accuracy based on a variety of different sized inputs.
>
> I do not think, that the effect of loop unrolling to the accuracy can be determined in general without considering the processed data.

True. However, my goals with MTIMESX were not necessarily to obtain a more accurate result, but to see if I could get a speed increase without a *decrease* in accuracy when compared to MATLAB. I found that loop unrolling was sometimes an effective way to do this (depends highly on system & compiler used).

> > I never saw the BLAS dot product result to be so bad that I would call it inaccurate.
>
> A short example, Matlab 2009a, Windows:
> x = ones(1, 3);
> y = [1e9, 1e-9, -1e9];
> x * y' % Calls BLAS
> >> 1e-9
> sum(x .* y)
> >> 0
>
> Whenever the BLAS dot product seems to be not accurate enough, try the XBLAS functions, which are stable and well tested also.

I think your earlier comment applies here also. If whatever you are doing with the data is sensitive enough that a method of summing is significantly affecting the results, you might question how the original data was processed to see if what you are summing even has the necessary resolution. As a side note, I think you need to be careful with the sum function since MATLAB has changed the method used for this somewhat recently (Bruno probably has more details on this).

(Maybe I will put different dot product method options, geared more towards accuracy improvement rather than speed improvement, in MTIMESX as a future enhancement.)

James Tursa

Subject: Double precision inner product

From: Derek O'Connor

Date: 20 Jan, 2011 01:22:06

Message: 19 of 22

"James Tursa" wrote in message <ih7m14$6ah$1@fred.mathworks.com>...
> "Derek O'Connor" wrote in message <ih71md$2e0$1@fred.mathworks.com>...
 --- snip --
> > This confirms what you say about Matlab's dot(x,y) -- slow and
> > innacurate (rel err = 1.6e-4)
> >
> > x'*y is fast but inaccurate (rel err = 3.4e-6)

> >
> > Derek O'Connor.
>
> I think you are comparing apples to oranges for some of these. In particular, labeling the x'*y calculation as inaccurate and the ddsumxy calculation as accurate. The x'y calculation may be doing intermediate calculations/storage in single precision whereas the ddsumxy method is doing all intermediate calculations/storage in double precision. In fact, the x'*y result is pretty close to being accurate to the full single precision mantissa length which is quite acceptable performance. So I don't see any accuracy issue with the x'*y method.
>
> I did extensive testing of dot product accuracy calculations when I coded my MTIMESX FEX submission, comparing the results to 100 decimal digit accurate calculations using the same inputs. I always saw the BLAS dot product routine to be among the best for accuracy, although as I mentioned previously you can sometimes get a slightly more accurate result with custom unrolled loops. The amount of unrolling (10) I eventually decided on was a compromise between speed and accuracy based on a variety of different sized inputs. I never saw the BLAS dot product result to be so bad that I would call it inaccurate.
>
> James Tursa

---------------------

James,

You are right about the accuracy of x'*y : it is as good as can be expected. I forgot it was doing the sum in single precision.

However, dot(x,y) seems to have lost 2 digits precision somewhere. It can't be the data because (x, y) = rand(10^6,1) is well-conditioned with respect to dot product. Higham says that trouble can be expected if abs(x'*y) << abs(x)'*abs(y). This condition does not hold in this case because x and y are positive.

Finally, the result of x(i)*y(i) must be twice the precision of x(i), y(i). Can this not be accumulated in a double precision sumxy which is rounded to single precision at the end of the loop? I thought this was always done, at least on IEEE FP machines.

Regards,

Derek

Subject: Double precision inner product

From: James Tursa

Date: 20 Jan, 2011 02:16:20

Message: 20 of 22

"Derek O'Connor" wrote in message <ih82nu$opr$1@fred.mathworks.com>...
>
> However, dot(x,y) seems to have lost 2 digits precision somewhere. It can't be the data because (x, y) = rand(10^6,1) is well-conditioned with respect to dot product. Higham says that trouble can be expected if abs(x'*y) << abs(x)'*abs(y). This condition does not hold in this case because x and y are positive.

I am not sure what the deal is with the dot function. I just recalled that it was slower and gave a different answer than a straightforward matrix multiply, and my somewhat limited investigations indicated that it was less accurate. So I tend to avoid it. In the m-file code for it the basic calculation is:

c = sum(conj(a).*b);

So it takes the explicit conjugate prior to the multiply, then allocates memory to hold the element-wise product, then calculates the sum. I can see how this takes more time than a simple matrix multiply. And the accuracy is going to depend on how the sum function works, and of course the fact that the intermediate results for the conj(a).*b are single precision. Other methods such as a BLAS routine may hold some of the intermediate results in higher precision and do some block summing. Bottom line is I don't really see why I would want to use the dot function for simple dot product calculations. The only advantage it seems is the ability to do the dot product along different dimensions, which I have never needed.

> Finally, the result of x(i)*y(i) must be twice the precision of x(i), y(i). Can this not be accumulated in a double precision sumxy which is rounded to single precision at the end of the loop? I thought this was always done, at least on IEEE FP machines.

I don't know. I think this depends on the processor and compiler being used. Certainly one can *force* double precision calculation and accumulation by explicitly coding it (ala Bruno's earlier post), but whether the compiler will, on its own, turn a single precision calculation/accumulation into a higher precision calculation/accumulation (e.g., in a register) is probably hard to control unless there are specific compiler settings available for this.

James Tursa

Subject: Double precision inner product

From: James Tursa

Date: 20 Jan, 2011 20:18:04

Message: 21 of 22

"James Tursa" wrote in message <ih7t5o$mup$1@fred.mathworks.com>...
> "Jan Simon" wrote in message <ih7nc4$4ma$1@fred.mathworks.com>...
> >
> > A short example, Matlab 2009a, Windows:
> > x = ones(1, 3);
> > y = [1e9, 1e-9, -1e9];
> > x * y' % Calls BLAS
> > >> 1e-9
> > sum(x .* y)
> > >> 0
> >
> > Whenever the BLAS dot product seems to be not accurate enough, try the XBLAS functions, which are stable and well tested also.
>
> I think your earlier comment applies here also. If whatever you are doing with the data is sensitive enough that a method of summing is significantly affecting the results, you might question how the original data was processed to see if what you are summing even has the necessary resolution. As a side note, I think you need to be careful with the sum function since MATLAB has changed the method used for this somewhat recently (Bruno probably has more details on this).
>
> James Tursa

Follow-up. It appears the change to sum and mtimes was made in version R2010a. Using 32-bit WinXP on an Intel Core 2 Duo machine:

L = 1e9;
s = 1e-9;
y = [1;1;1];

x = [L;s;-L];
R2006b - R2009b
x'*y --> 1e-9
sum(x.*y) --> 0
R2010a - R2010b
x'*y --> 0
sum(x.*y) --> 0

x = [L;-L;s];
R2006b - R2009b
x'*y --> 0
sum(x.*y) --> 1e-9
R2010a - R2010b
x'*y --> 1e-9
sum(x.*y) --> 1e-9

x = [s;L;-L];
R2006b - R2009b
x'*y --> 0
sum(x.*y) --> 0
R2010a - R2010b
x'*y --> 0
sum(x.*y) --> 0

And it would not surprise me to find that these results could change from processor to processor as well. So you can't rely on mtimes and/or sum to do these types of calculations to the "best" precision. That being said, I would reiterate my previous statements that if you need this type of calculation to give you a 1e-9 result then I would seriously question the original data. Is that 1e9 really exactly 1e9 from some calculation where you know the trailing bits beyond the actual stored mantissa are really exactly zero? If not, and this is some real world data that is being summed where the 1e9 is some type of measurement, then both results shown above (0 and 1e-9) are really just garbage because of the massive cancellation involved in the calculation. So while the 1e-9 answer certainly is pleasing to the eye, it is no better than the 0 answer unless you know some very specific things
about the input data. Again, I am not too worried about the accuracy of mtimes in the typical cases I use it for.

James Tursa

Subject: Double precision inner product

From: Jan Simon

Date: 21 Jan, 2011 10:07:04

Message: 22 of 22

Dear James,

> And it would not surprise me to find that these results could change from processor to processor as well. So you can't rely on mtimes and/or sum to do these types of calculations to the "best" precision.

I expect calls to the BLAS functions to be "best" applicable method, as long as XBLAS is too slow.
The [1e9, 1e-9, -1e9] data are artificial. If a "real world problem" contains such data, the sum is simply instable. A scientific colution of this problem must include a "meaningful" variation of the inputs and the instability will be detected in a very high sensitivity of the outputs to a tiny variation of the input.

I've processed a real world problems, which is affected by the limited accuracy of the dot product: The estimation of the Hessian for solving an ODE with internal numerical differentiation. The sensitivity of the simulated physical system appears as massive cancellation error in the centered differences to estimate the 2nd derivatives.
But it is always trivial to handle such problems: The result is not just the final number, but it must be accompanied by the Wronski-matrix (sensitivity of outputs to variation of inputs). E.g. "the trajectory of a pencil standing on the tip" can be calculated easily, but without adding "but the faintest disturbance will lead to a completely different answer" any solution is hillarious.

However, it can be a hard problem already to find a "meaningful" variation of the inputs! For the pencil the measurement error of the initial position and velocity is fine, and the rotation of the earth can be ignored. For a system of GPS satellites it is much more complicated...
For "real world problems" the limited accuracy of the inner product can be controlled easily (but not avoided), but the limited accuracy of the parameterization of the phsyical (chemical, biological, etc) problem remains a combination of science, intuition, wild guessing and intented manipulation of the results.

Kind regards, Jan

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us