Outer prod
Marcus M. Edvall
>> a = [1 2 3 4 5]';
>> b = a/1e-5;
>> a*b'-(a*b')'
ans =
1.0e-009 *
0 0 0.0582 0 0
0 0 0.1164 0 0
-0.0582 -0.1164 0 -0.2328 -0.2328
0 0 0.2328 0 0
0 0 0.2328 0 0
Best wishes, Marcus
Tomlab Optimization Inc.
http://tomopt.com
Loren Shure
edv...@gmail.com says...
Answer from Cleve:
a*b' is symmetric if and only if b is an exact scalar multiple of a. In
this case, b is not an exact scalar multiple of a.
Tim Davis
ans =
1.0e-016 *
0 0 0.1058 0 0
0 0 0.2117 0 0
-0.1058 -0.2117 0 -0.4233 -0.4233
0 0 0.4233 0 0
0 0 0.4233 0 0
The relative error is O(eps) ... it's just roundoff error.
No surprise here.
Michael Hosea
might argue that there should be no roundoff error here because everything
is a floating point integer, and all the floating point multiplications and
subtractions should be done exactly with this data. The problem here is
caused by b = a/1e-5. Note that
>> 1/1e-5 == 1e5
ans =
0
Because of this, b does not end up being an array of floating point
integers. However, if we trim off the noise first:
>> b = round(b)
b =
100000
200000
300000
400000
500000
>> a*b' - (a*b')'
ans =
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
--
Mike
Tim Davis
<fcalm4$8r6$1...@fred.mathworks.com>...
Mike,
That assumes that one should calculate x/y as x*(1/y). In
this case, you are right; 1/1e-5 is more accurately computed
as 1*1e5. However, the number 1e-5 itself is not
representable as a floating-point number. The MATLAB
interpreter would have to symbolically determine that x/y
could be done as x*(1/y) where now 1/y *is* a representable
floating-point number (and in particular a "flint", or an
integer stored as a floating-point number).
That's a tall order, and I don't think it's worth it.
MATLAB uses floating-point, except in the Symbolic Toolbox,
which means it wants speed and accepts a little
floating-point roundoff as a tradeoff.
x*(1/y) is multiplying by the inverse of y, which is
*always* a bad idea except for peculiar cases like this when
1/y is already known *and* known to be more accurate. In
general x*(1/y) results in twice the work and twice the
roundoff error. x*(1/y) if y is denormal would give Inf
instead of a correct result.
If a MATLAB user wants to do b/1e-5, then b*1e5 is better
... so just code it that way (assuming 1e5 is a hard
constant), rather than use round(b). If on the other hand
you are doing x/y where y could be anything, then don't try
x*(1/y).
Michael Hosea
I don't understand why you wrote a reply to me as if to argue some point. I
was merely pointing out the source of the actual numerical error in this
case, not arguing that it shouldn't be there. I did not think it was
obvious where the error was really coming from.
--
Mike
"Tim Davis" <da...@cise.ufl.edu> wrote in message
news:fcb7g2$843$1...@fred.mathworks.com...
Marcus M. Edvall
is symmetric, i.e. something else than a*b'?
It makes no sense to me why
a(1)*b(3) ~= b(3)*a(1)
Since it was symmetric before Matlab 7.2 it would be interesting to
understand why the decision has been made that it should not
necessarily be symmetric.
I'd hate to write a C MEX to multiply two vectors.
c=1/2*(a*b'+(a*b')') would be a bit silly. Hmmm.
Michael Hosea
symmetric (e.g.
>> [1;1]*[-1,1]
ans =
-1 1
-1 1
For a*b' to be symmetric requires that
a(i)*b(j) == b(i)*a(j),
which in your example is a(1)*b(3) == a(3)*b(1). One might expect this to
be satisfied because b = a/1e-5. The requirement for symmetry works out to
be
a(1)*(a(3)/1e-5) == a(3)*(a(1)/1e-5)
But you don't generally get the same floating point results down to the last
bit when you regroup floating point operations. In particular,
>> format hex
>> 1*(3/1e-5)
ans =
41124f8000000000
>> 3*(1/1e-5)
ans =
41124f7fffffffff
We should expect to get errors on the order of eps(a(i)*b(j)). If it
happens to work out nicely in some cases, e.g.
>> 2*(3/0.1)
ans =
404e000000000000
>> 3*(2/0.1)
ans =
404e000000000000
so much the better.
--
Mike
"Marcus M. Edvall" <edv...@gmail.com> wrote in message
news:1189704805.9...@w3g2000hsg.googlegroups.com...
Michael Hosea
the behavior of the software hasn't changed in this regard. Perhaps in an
earlier version of MATLAB you used a scaling factor that happened to work,
e.g.
>> b = a/1e-4
b =
10000
20000
30000
40000
50000
>> a*b' - (a*b')'
ans =
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
As far as forcing exact symmetry goes, A = (A+A.')*0.5 is something I've
been known to do. I don't know a faster way.
--
Mike
"Michael Hosea" <mho...@mathworks.com> wrote in message
news:fcc0kr$o6u$1...@fred.mathworks.com...
Tim Davis
> Tim,
>
> I don't understand why you wrote a reply to me as if to
argue some point. I
> was merely pointing out the source of the actual numerical
error in this
> case, not arguing that it shouldn't be there. I did not
think it was
> obvious where the error was really coming from.
> --
> Mike
>
Sorry ... I didn't mean to argue. I thought you raised a
great point, which is why I referred to your post. I was
just pointing out that it would be difficult for MATLAB, by
itself, to fix the problem along those lines. It would
require the recoding of b=b/1e-5 to b=b*1e5, by the user.
I'm sorry if I sounded like I was arguing.
Thanks,
Tim
Michael Hosea
I was in a hurry to go off to lunch at the time. :)
I completely agree. It's fundamental, a predictable (if not always
predicted) consequence of using binary floating point. The confusing thing
about this particular problem is that it's easy to see a/1e-5 and think this
is the same as a*1e5, but of course it is not.
--
Mike
"Tim Davis" <da...@cise.ufl.edu> wrote in message
news:fccdi3$9k2$1...@fred.mathworks.com...
Tim Davis
That's the dilemma of a newsgroup posting, and email too ...
it can sound so formal. Without, say, poetry ;-), it's hard
to express emotion or nuance.
... don't worry, I'm not about to break into rhyme (although
"flint" would make a great word to use in a poem or pun!) ...
I would guess that (A+A')/2 is the best technique to ensure
symmetry (or equally (A+A')*.5, since .5 truly is 1/2 in
floating point...). Alternatively, one could do A =
tril(A)+tril(A,-1)'. Not sure which is faster; the latter
at least does no flops. For the complex case ... I suppose
it depends on whether you want symmetry or Hermitianity (is
that a word ??).
Bobby Cheng
The fastest way I can think of is
(a*a')./1e-5
MATLAB ensures that a*a' returns Hermitian result for a few release now.
---Bob.
"Tim Davis" <da...@cise.ufl.edu> wrote in message
news:fccovl$45$1...@fred.mathworks.com...
Peter Perkins
> Is there any way after Matlab 7.2 to calculate an outer product that
> is symmetric, i.e. something else than a*b'?
>
> It makes no sense to me why
>
> a(1)*b(3) ~= b(3)*a(1)
>
> Since it was symmetric before Matlab 7.2 it would be interesting to
> understand why the decision has been made that it should not
> necessarily be symmetric.
>
> I'd hate to write a C MEX to multiply two vectors.
>
> c=1/2*(a*b'+(a*b')') would be a bit silly. Hmmm.
Just to be clear: As Bobby points out, a*a' (and indeed A*A' and A'*A) is
computed in such a way as to be exactly symmetric, has been for several
releases. So if what you really want is alpha*(a*a'), then it is not necessary
to write c=1/2*(a*a'+(a*a')'), because a*a' _will be_ symmetric. As Mike
describes, the problem in your example was not the outer product, but rather the
original calculation of b, or more to the point, the regrouping of floating
point operations.
Hope this helps.