Some unstable phenomenon of matrix operation, what is the reason?
Sui Huang
numerical error, they should make the identical error...), but during
the computation they are not. I'm confused, where is the difference
from? Can anybody tell me?
Thank you!
Sui
[code]
K>> H=Pts(:,T)'*Pts(:,T);
K>> H
H =
1.0000 0.1842 -0.2822 -0.3202 0.4244 -0.2148
-0.5205 -0.2885 0.2672
0.1842 1.0000 -0.4791 -0.3465 -0.3734 0.0292
-0.6650 -0.2129 -0.4371
-0.2822 -0.4791 1.0000 0.5782 0.3820 -0.1440
0.2836 -0.2338 0.1222
-0.3202 -0.3465 0.5782 1.0000 0.4099 0.3655
-0.0957 -0.4549 0.5104
0.4244 -0.3734 0.3820 0.4099 1.0000 0.0152
-0.1602 -0.2325 0.8239
-0.2148 0.0292 -0.1440 0.3655 0.0152 1.0000
-0.2267 0.1841 0.3717
-0.5205 -0.6650 0.2836 -0.0957 -0.1602 -0.2267
1.0000 0.5893 -0.2490
-0.2885 -0.2129 -0.2338 -0.4549 -0.2325 0.1841
0.5893 1.0000 -0.1927
0.2672 -0.4371 0.1222 0.5104 0.8239 0.3717
-0.2490 -0.1927 1.0000
K>>
K>> testA=Pts(:,T);
K>> testH=testA'*testA;
K>>
K>> testH-H
ans =
1.0e-015 *
0 0 0 0 0 0
0 0 0
0 0 0.1110 0.1110 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 0 0 0 0
0 -0.0555 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
0 0 0
0 0 0 0 0 0
0 0 0
K>>
[\code]
Christopher Hulbert
I would venture to guess that 1e-16 is numerical noise and is less than eps which all the guaranteed
accuracy you can expect in floating point arithmetic.
Sui
> than eps which all the guaranteed
> accuracy you can expect in floating point arithmetic.
Yes, for the posted case the noise can be ignored. But if the
size of matrix is increased to larger size such as 100 or 1000, then
the noise will be significant.
In my experiment, I used a H of dimension 10 as the argument of
quadratic optimization function in optimization toolbox. Then it
returns warning that H is not symmetric.
Though I can use some quick way to solve this problem, but I
could not understand where is the noise from...
Sui
Peter Boettcher
> H and testH should be exactly the same by intuition (even if there is
> numerical error, they should make the identical error...), but during
> the computation they are not. I'm confused, where is the difference
> from? Can anybody tell me?
> Thank you!
>
> Sui
>
> [code]
> K>> H=Pts(:,T)'*Pts(:,T);
> K>> H
>
> H =
>
[snip]
> K>>
> K>> testA=Pts(:,T);
> K>> testH=testA'*testA;
> K>>
> K>> testH-H
>
> ans =
>
> 1.0e-015 *
[snip a few small entries]
The only thing I can come up with is that MATLAB automatically detects
A'*A operations, and computes a symmetric output. Perhaps that
detection is somehow disrupted by the subscripting, so the full output
is computed.
Or the accelerator does something funny.
Or order of computation means that more numbers stay in the FPU
extended precision registers (80-bit on Intel) rather than being
stored back to 64-bit memory representation.
I'm only guessing here.
--
Peter Boettcher <boet...@ll.mit.edu>
MIT Lincoln Laboratory
MATLAB FAQ: http://www.mit.edu/~pwb/cssm/
Roger Stafford
<super...@hotmail.NOSPAM.com> wrote:
--------------------
Hello Sui,
I have not been able to duplicate your results with my matlab version
(4a), though I tried many tens of thousands of randomly determined A's and
T's. Can you tell a little more about T? I take it T is a nine-element
vector of indices. What is the number of columns in Pts? Are the indices
of T in ascending order? Do they repeat? Are you absolutely sure they
are all exact integers in the proper range?
If there is any difference in the sequencing of operations between the
two matrix multiplications, Pts(:,T)'*Pts(:,T) and testA'*testA, that
could introduce differences in round-off errors. (For example, with three
floating point scalars a, b, and c, a*(b*c) and (a*b)*c can easily yield
different results, though only out at the least bit.) However, it is
difficult for me to imagine how that could apply to your case. When
matlab is presented with an expression, Pts(:,T), it is supposed to first
calculate its value, just as it did in the assignment testA=Pts(:,T), and
only then proceed with the matrix multiplication, and that ought then to
produce identical results to using the precalculated testA. To do
otherwise would suggest some kind of bug in MathWorks' matrix
multiplication logic in your version.
If you are very sure of your results, I would suggest sending them to
MathWorks' technical support for their scrutiny.
(Remove "xyzzy" and ".invalid" to send me email.)
Roger Stafford
John D'Errico
<ellieandrogerxyzz...@pool0726.cvx22-bradley.dialup.earth
link.net>,
ellieandr...@mindspring.com.invalid (Roger Stafford) wrote:
> I have not been able to duplicate your results with my matlab version
> (4a), though I tried many tens of thousands of randomly determined A's and
> T's. Can you tell a little more about T? I take it T is a nine-element
> vector of indices. What is the number of columns in Pts? Are the indices
> of T in ascending order? Do they repeat? Are you absolutely sure they
> are all exact integers in the proper range?
A release or so ago they added a speed enhancement
where the parser recognizes A'*A operations. It takes
advantage of the special form, not needing to generate
A', so there is a speed & memory savings there, plus
probably a speed gain by only doing 1/2 as much work.
Odds are it does not pick up that this product is of
the same form. I'll guess that
feature accel off
might confirm this fact.
Simplest is to force symmetry on the result if there
is a problem.
John
--
The best material model of a cat is another, or
preferably the same, cat.
A. Rosenblueth, Philosophy of Science, 1945
Sui Huang
Thank you for your hard testing. Sorry about the late reply.
> I have not been able to duplicate your results with my matlab
> version
> (4a), though I tried many tens of thousands of randomly determined
A's and T's.
I'm using Matlab 7 SP2, this may be a reason of difference.
> Can you tell a little more about T? I take it T is a
> nine-element vector of indices. What is the number of
> columns in Pts?
Yes, T is a row vector of indices, though the type is double, but
still have exactly integer values, which can be use as column
indices. If length of T is d, then number of columns in Pts is d^2.
The number of columns in Pts is very large. Also if d increases, the
pheomenone will be more obvious. d=9 is the smallest size that I can
find the pheomenone to post it easily.
>Are the indices of T in ascending order?
Yes
> Do they repeat?
No
> Are you absolutely sure they are all exact integers in the proper
range?
Yes, I use T as indices for a long time (during a year, I wrote many
routines, indexing matrix columns with T, and ran a lot of tests for
something else.)
Thank you for your suggestion!
Sui
Sui Huang
> where the parser recognizes A'*A operations. It takes
> advantage of the special form, not needing to generate
> A', so there is a speed & memory savings there, plus
> probably a speed gain by only doing 1/2 as much work.
>
> Odds are it does not pick up that this product is of
> the same form. I'll guess that
>
> feature accel off
>
> might confirm this fact.
>
> Simplest is to force symmetry on the result if there
> is a problem.
John:
As you said, the parser maybe the reason. I just tried a simpler
example. Seems MATLAB deals with the A(index_lists)'*A(index_list) in
a special way. However, I dont understand why they do this, becasue
when I experiment it in large size (1000 by 1000) in the example, it
even slows down the routine. This may be a bug, I guess...
Sui
[code]
>> A=rand(8)
A =
0.5828 0.5298 0.7942 0.7680 0.3200 0.6833
0.3705 0.6831
0.4235 0.6405 0.0592 0.9708 0.9601 0.2126
0.5751 0.0928
0.5155 0.2091 0.6029 0.9901 0.7266 0.8392
0.4514 0.0353
0.3340 0.3798 0.0503 0.7889 0.4120 0.6288
0.0439 0.6124
0.4329 0.7833 0.4154 0.4387 0.7446 0.1338
0.0272 0.6085
0.2259 0.6808 0.3050 0.4983 0.2679 0.2071
0.3127 0.0158
0.5798 0.4611 0.8744 0.2140 0.4399 0.6072
0.0129 0.0164
0.7604 0.5678 0.0150 0.6435 0.9334 0.6299
0.3840 0.1901
>>
>> H=A'*A;
>> H-H'
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 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 0 0 0 0 0 0
0 0 0 0 0 0 0 0
>>
>> testH=A(:,:)'*A(:,:);
>> testH'-testH
ans =
1.0e-015 *
0 0 0.2220 0 0 0
0 0
0 0 0 0 0 0
0 0
-0.2220 0 0 0 0 -0.2220
0 0
0 0 0 0 0 0
0 0
0 0 0 0 0 0
0 0
0 0 0.2220 0 0 0
0 0
0 0 0 0 0 0
0 0
0 0 0 0 0 0
0 0
>>
>> A=rand(1000);
>> tic; H=A'*A; toc;
Elapsed time is 0.701011 seconds.
>> tic; H=A'*A; toc;
Elapsed time is 0.755078 seconds.
>> tic; H=A'*A; toc;
Elapsed time is 0.753196 seconds.
>>
>>
>> tic; testH=A(:,:)'*A(:,:); toc;
Elapsed time is 1.317488 seconds.
>> tic; testH=A(:,:)'*A(:,:); toc;
Elapsed time is 1.270066 seconds.
>> tic; testH=A(:,:)'*A(:,:); toc;
Elapsed time is 1.277384 seconds.
>>
>>
>> tic; testH=A(:,1:1000)'*A(:,1:1000); toc;
Elapsed time is 1.261370 seconds.
>> tic; testH=A(:,1:1000)'*A(:,1:1000); toc;
Elapsed time is 1.249645 seconds.
>> tic; testH=A(:,1:1000)'*A(:,1:1000); toc;
Elapsed time is 1.195533 seconds.
>>
>>
>> tic; testA=A(:,:); testH=testA'*testA; toc;
Elapsed time is 0.769432 seconds.
>> tic; testA=A(:,:); testH=testA'*testA; toc;
Elapsed time is 0.712737 seconds.
>> tic; testA=A(:,:); testH=testA'*testA; toc;
Elapsed time is 0.758192 seconds.
>>
>>
>>
>> tic; testA=A(:,:); tran_testA=A(:,:)'; testH=testA'*testA; toc;
Elapsed time is 0.847046 seconds.
>> tic; testA=A(:,:); tran_testA=A(:,:)'; testH=testA'*testA; toc;
Elapsed time is 0.832330 seconds.
>> tic; testA=A(:,:); tran_testA=A(:,:)'; testH=testA'*testA; toc;
Elapsed time is 0.825600 seconds.
>>
[/code]
Bobby Cheng
can recogise H=Pts(:,T)'*Pts(:,T); is in fact A = Pts(:,T); H = A*A';. and
optimize this behind the scene. This is for the future and I wish MATLAB
could do that too.
For now, A*A' is the easiest case to detect and is very useful to ensure the
product is Hermitian.
The following code pattern is common and can be simplified using this fact
from
Y = A*A'; Y = (Y+Y')/2; to Y = A*A';
---Bob.
"Sui Huang" <super...@hotmail.NOSPAM.com> wrote in message
news:ef12...@webx.raydaftYaTP...
Peter Boettcher
> One can say that this is a limitation of MATLAB parser. Maybe one day MATLAB
> can recogise H=Pts(:,T)'*Pts(:,T); is in fact A = Pts(:,T); H = A*A';. and
> optimize this behind the scene. This is for the future and I wish MATLAB
> could do that too.
>
> For now, A*A' is the easiest case to detect and is very useful to ensure the
> product is Hermitian.
>
> The following code pattern is common and can be simplified using this fact
> from
>
> Y = A*A'; Y = (Y+Y')/2; to Y = A*A';
I've mumbled about this sort of thing before. Maybe you can explain
it to me:
I've often wondered why MATLAB does not have delayed-operation
transposes. That is, A' would simply create a linked copy of A with a
'transpose' flag set. The old API to access the data itself (mxGetPr,
or whatever you use internally) would create the actual transpose, so
that everything would still work without recoding. Then a new API
would optionally retain that flag
(mxGetPrEx(... SUPPORT_FLAGS_TRANSPOSE)) so that the low-level stuff
could gradually support it. BLAS operations natively support
transpose, so that's an obvious place to start. A*B' could be done
without ever creating a temporary B'.
This makes it really easy to recognize A*A', and add a 'symmetric'
flag, so that solvers wouldn't have to 'detect' symmetric
matrices...
The possibilities are endless... banded storage, etc.
Sui Huang
**********************
In your code below, the reason why line 2 takes longer than line 3 is
because in
line 2 there is a function call to a MATLAB built-in function,
SUBSREF. SUBSREF
is designed to be used by the MATLAB interpreter to handle indexed
references to
objects. It is an overloaded method for A(I), A{I} and A.field, where
A is an
array.
1. A=rand(1000);
2. tic; H=A(:,:)'*A(:,:); toc;
3. tic; tA=A; tH=tA'*tA; toc;
4. any(any(tH'-tH))
5. any(any(H'-H))
Line 4 and line 5 resulted in the correct answer: zero. I was unable
to
reproduce the error that you are encountering.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Is there anybody can reproduce the error?
Thank you!
Bobby Cheng
transpose.
The problem here is Pts(:,T)' and Pts(:,T) actually created two different
copy of the same matrix.
So we are talking about a special case of A*B' where A and B are two copys
of the same things. The only way to know the two matrices are actually
equal, one has to check every single element. This will slow down the
general case in which A and B are different.
A nice thing to do here is for MATLAB to recognise that the user is trying
to create the same things twice and be smart about this. This is an extra
level of complexity for the parser that MATLAB don't currently have.
Does this make things clearer? I hope I have not confused the issue.
---Bob.
"Peter Boettcher" <boet...@ll.mit.edu> wrote in message
news:m38xyga...@coyote.llan.ll.mit.edu...
Peter Boettcher
[Reordering top-post]
"Bobby Cheng" <bch...@mathworks.com> writes:
>
> "Peter Boettcher" <boet...@ll.mit.edu> wrote in message
> news:m38xyga...@coyote.llan.ll.mit.edu...
[snip]
>>
>> I've often wondered why MATLAB does not have delayed-operation
>> transposes. That is, A' would simply create a linked copy of A with a
>> 'transpose' flag set. The old API to access the data itself (mxGetPr,
>> or whatever you use internally) would create the actual transpose, so
>> that everything would still work without recoding. Then a new API
>> would optionally retain that flag
>> (mxGetPrEx(... SUPPORT_FLAGS_TRANSPOSE)) so that the low-level stuff
>> could gradually support it. BLAS operations natively support
>> transpose, so that's an obvious place to start. A*B' could be done
>> without ever creating a temporary B'.
>>
>> This makes it really easy to recognize A*A', and add a 'symmetric'
>> flag, so that solvers wouldn't have to 'detect' symmetric
>> matrices...
>>
>> The possibilities are endless... banded storage, etc.
>
>
> We have already done this sort of things. A*B' does not actually perform the
> transpose.
>
> The problem here is Pts(:,T)' and Pts(:,T) actually created two different
> copy of the same matrix.
>
> So we are talking about a special case of A*B' where A and B are two copys
> of the same things. The only way to know the two matrices are actually
> equal, one has to check every single element. This will slow down the
> general case in which A and B are different.
>
> A nice thing to do here is for MATLAB to recognise that the user is trying
> to create the same things twice and be smart about this. This is an extra
> level of complexity for the parser that MATLAB don't currently have.
>
> Does this make things clearer? I hope I have not confused the issue.
That's not my point. I'm suggesting to do it at a level that does not
rely on the parser detecting certain text patterns. The case of
Pts(:,T) could still not be solved without additional smarts. The way
I describe would handle such things as:
B=A';
A*B;
as a symmetric product with zero overhead. And would allow
zero-overhead simplification of stuff like sum(A.'), because "sum"
could see the transpose flag, and simply execute sum(A,2).' And would
know that transpose of a matrix with the SYMMETRIC flag set is simply
itself.
Pts(:,T) would have to be done at a higher level, recognizing the
common expression Pts(:,T). Once that's done, you're all set.
temp1=Pts(:,T);
Then temp1'*temp1 is evaluated normally, with the transpose flag.
Compiler optimizers have been extracting common expressions for a very
long time. I'm sure the MATLAB accelerator can support that sort of
optimization, or should.
Whatever. The parser approach is helping some. I'm not holding my
breath for this sort of improvement.