Some unstable phenomenon of matrix operation, what is the reason?

Discussion in 'MATLAB' started by Christopher Hulbert, Aug 29, 2005.

1. Christopher HulbertGuest

Christopher Hulbert, Aug 29, 2005

2. Sui HuangGuest

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 =

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]```

Sui Huang, Aug 29, 2005

3. SuiGuest

I would venture to guess that 1e-16 is numerical noise and is less
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

Sui, Aug 29, 2005
4. Peter BoettcherGuest

Peter Boettcher, Aug 29, 2005
5. Roger StaffordGuest

Roger Stafford, Aug 29, 2005
6. John D'ErricoGuest

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

John D'Errico, Aug 29, 2005
7. Sui HuangGuest

Hi! Roger:
A's and T's.
I'm using Matlab 7 SP2, this may be a reason of difference.
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.
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.)

Sui

Sui Huang, Sep 1, 2005
8. Sui HuangGuest

A release or so ago they added a speed enhancement
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 =

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

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

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
```

Sui Huang, Sep 1, 2005
9. Bobby ChengGuest

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';

---Bob.

Bobby Cheng, Sep 1, 2005
10. Peter BoettcherGuest

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.

Peter Boettcher, Sep 1, 2005
11. Sui HuangGuest

Here is the reply of technical support:

**********************
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!

Sui Huang, Sep 1, 2005
12. Bobby ChengGuest

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
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.

Bobby Cheng, Sep 1, 2005
13. Peter BoettcherGuest

[Reordering top-post]

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.

Peter Boettcher, Sep 2, 2005