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:
a.*a ~= a*a'

Subject: a.*a ~= a*a'

From: Michael

Date: 2 Apr, 2009 01:55:06

Message: 1 of 8

Here's something interesting.

>> p = [-1:-.1:-100];
>> sum(p .* p) - sum(p * p')
ans =
     4.656612873077393e-10
>> eps(sum(p*p'))
ans =
     4.656612873077393e-10

The sum of element-by-element multiplication is different than matrix multiplication. It's only off here by eps, but in my tests for other kinds of values I get larger deviations.

Seems it's probably because of the algorithms working in different order and thereby having different float rounding error accumulation:

>> sum(p .* p) - sum(p(991:-1:1) * p(991:-1:1)')
ans =
     0

Interesting. Anyone know for sure if this is the reason?

Cheers,
Michael

Subject: a.*a ~= a*a'

From: Matt Fig

Date: 2 Apr, 2009 03:12:00

Message: 2 of 8

See this thread.
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245808

Subject: a.*a ~= a*a'

From: Roger Stafford

Date: 2 Apr, 2009 03:23:01

Message: 3 of 8

"Michael" <michael@circular-logic.com> wrote in message <gr15tq$f9b$1@fred.mathworks.com>...
> Here's something interesting.
>
> >> p = [-1:-.1:-100];
> >> sum(p .* p) - sum(p * p')
> ans =
> 4.656612873077393e-10
> >> eps(sum(p*p'))
> ans =
> 4.656612873077393e-10
>
> The sum of element-by-element multiplication is different than matrix multiplication. It's only off here by eps, but in my tests for other kinds of values I get larger deviations.
>
> Seems it's probably because of the algorithms working in different order and thereby having different float rounding error accumulation:
>
> >> sum(p .* p) - sum(p(991:-1:1) * p(991:-1:1)')
> ans =
> 0
>
> Interesting. Anyone know for sure if this is the reason?
>
> Cheers,
> Michael

  The numbers which matlab uses are binary floating point numbers which are incapable of representing -.1 exactly. This does not in itself cause the errors you are seeing but it does cause the numbers in your 'p' to have strings of both 1's and 0's running all the way down to the least bit position. This in turn can lead to round-off errors with such numbers that can produce slightly different results for what would ordinarily be regarded as equivalent operations.

  Here is a problem that you can work by hand to demonstrate this to your own satisfaction. Do two different multiplications in four-digit binary numbers with a = .1100, b = .1101, and c = .1110 where these are binary fractions, not decimal, and always round off your answer to the nearest four digits at each step.

 (a*b)*c will produce .1001 and
 a*(b*c) will give .1000

Of course neither answer is exactly correct, but the associative law of multiplication has thereby been violated. The same thing can happen with the 53-bit binary significands matlab uses.

Roger Stafford

Subject: a.*a ~= a*a'

From: James Tursa

Date: 2 Apr, 2009 15:38:01

Message: 4 of 8

"Michael" <michael@circular-logic.com> wrote in message <gr15tq$f9b$1@fred.mathworks.com>...
> Here's something interesting.
>
> >> p = [-1:-.1:-100];
> >> sum(p .* p) - sum(p * p')
> ans =
> 4.656612873077393e-10
> >> eps(sum(p*p'))
> ans =
> 4.656612873077393e-10
>
> The sum of element-by-element multiplication is different than matrix multiplication. It's only off here by eps, but in my tests for other kinds of values I get larger deviations.
>
> Seems it's probably because of the algorithms working in different order and thereby having different float rounding error accumulation:
>
> >> sum(p .* p) - sum(p(991:-1:1) * p(991:-1:1)')
> ans =
> 0
>
> Interesting. Anyone know for sure if this is the reason?
>
> Cheers,
> Michael

I will add my 2 cents to this.

If you are on a PC, then there is a numeric co-processor involved that does calculations in 80-bit arithmetic. So the sum(p .* p) formulation does each individual element-by-element calculation by converting each operand to 80-bit, then doing the multiply, then converting each individual result back to 64-bit for storage in the intermediate vector result. So an intermediate storage vector the size of p is needed. Then the sum is done, again likely by converting to 80-bit, doing the sum, then converting back to 64-bit.

For the sum(p * p') formulation (you don't need the sum() function here, btw), the p' is formed as a simple reshape of p with no data copying. You can see this as follows:

format debug
a = rand(1,5)
b = a'

You can see that the pr of both variables is the same.

Then the p * p' calculation is done by calling a BLAS routine, which again probably converts to 80-bit in the background, but in this case the actual summing may be done in 80-bit as well. i.e., the individual element-by-element multiplies are kept in 80-bit as they are added to the sum. This gives you a chance for lower cancellation errors in the sum calculation (whether you actually get much benefit here would depend on the actual data).

In short, p * p' is the better formulation for this particular calculation because it 1) doesn't require a large intermediate vector result in the calculation, 2) it probably gives you a better chance for smaller cancellation errors in the result, at least on a PC, and 3) It is faster, largely because of 1). e.g.,

>> a = rand(1,10000000);
>> tic;sum(a.*a);toc
Elapsed time is 0.273004 seconds.
>> tic;a*a';toc
Elapsed time is 0.029888 seconds.

James Tursa

Subject: a.*a ~= a*a'

From: Michael

Date: 2 Apr, 2009 22:19:01

Message: 5 of 8

"James Tursa" <aclassyguywithaknotac@hotmail.com> wrote in message <gr2m4p$c5b$1@fred.mathworks.com>...

>
> If you are on a PC, then there is a numeric co-processor involved that does calculations in 80-bit arithmetic. So the sum(p .* p) formulation does each individual element-by-element calculation by converting each operand to 80-bit, then doing the multiply, then converting each individual result back to 64-bit for storage in the intermediate vector result. So an intermediate storage vector the size of p is needed. Then the sum is done, again likely by converting to 80-bit, doing the sum, then converting back to 64-bit.
>
> For the sum(p * p') formulation (you don't need the sum() function here, btw), the p' is formed as a simple reshape of p with no data copying. You can see this as follows:
>
> format debug
> a = rand(1,5)
> b = a'
>
> You can see that the pr of both variables is the same.
>
> Then the p * p' calculation is done by calling a BLAS routine, which again probably converts to 80-bit in the background, but in this case the actual summing may be done in 80-bit as well. i.e., the individual element-by-element multiplies are kept in 80-bit as they are added to the sum. This gives you a chance for lower cancellation errors in the sum calculation (whether you actually get much benefit here would depend on the actual data).
>

Thanks James, this makes sense. I hadn't thought about Matlab using BLAS in one instance and not in the other, which could be compiled differently even if the algorithm otherwise is similar. And the part about the intermediate vector possibly getting converted back to 64-bit before the summation is very helpful too. That makes sense as an additional source for different fp rounding errors. I'm on an Intel Mac so I'm assuming it has the 80-bit coprocessor.

Cheers,
Michael

Subject: a.*a ~= a*a'

From: Steve Amphlett

Date: 3 Apr, 2009 06:59:02

Message: 6 of 8

"James Tursa" <aclassyguywithaknotac@hotmail.com> wrote in message

> For the sum(p * p') formulation (you don't need the sum() function here, btw), the p' is formed as a simple reshape of p with no data copying. You can see this as follows:
>
> format debug
> a = rand(1,5)
> b = a'
>
> You can see that the pr of both variables is the same.

For a 1D column/row vector, yes. But not in a 2D array surely.

Subject: a.*a ~= a*a'

From: James Tursa

Date: 3 Apr, 2009 15:23:02

Message: 7 of 8

"Steve Amphlett" <Firstname.Lastname@Where-I-Work.com> wrote in message <gr4c3m$6pt$1@fred.mathworks.com>...
> "James Tursa" <aclassyguywithaknotac@hotmail.com> wrote in message
>
> > For the sum(p * p') formulation (you don't need the sum() function here, btw), the p' is formed as a simple reshape of p with no data copying. You can see this as follows:
> >
> > format debug
> > a = rand(1,5)
> > b = a'
> >
> > You can see that the pr of both variables is the same.
>
> For a 1D column/row vector, yes. But not in a 2D array surely.

Yes, a 2D array transpose would require a data copy, of course. But that was not the problem OP posted.

James Tursa

Subject: a.*a ~= a*a'

From: Steve Amphlett

Date: 3 Apr, 2009 15:30:17

Message: 8 of 8

"James Tursa" <aclassyguywithaknotac@hotmail.com> wrote in message <gr59km$7it$1@fred.mathworks.com>...
> "Steve Amphlett" <Firstname.Lastname@Where-I-Work.com> wrote in message <gr4c3m$6pt$1@fred.mathworks.com>...
> > "James Tursa" <aclassyguywithaknotac@hotmail.com> wrote in message
> >
> > > For the sum(p * p') formulation (you don't need the sum() function here, btw), the p' is formed as a simple reshape of p with no data copying. You can see this as follows:
> > >
> > > format debug
> > > a = rand(1,5)
> > > b = a'
> > >
> > > You can see that the pr of both variables is the same.
> >
> > For a 1D column/row vector, yes. But not in a 2D array surely.
>
> Yes, a 2D array transpose would require a data copy, of course. But that was not the problem OP posted.

Transposing a 2D array in a MEX is a real PITA. Simple to do; simple to get wrong. Always wrong the first time.

Tags for 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