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

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

  1.  
    Christopher Hulbert, Aug 29, 2005
    #1
    1. Advertisements

  2. Christopher Hulbert

    Sui Huang Guest

    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
    #2
    1. Advertisements

  3. Christopher Hulbert

    Sui Guest

    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
    #3
  4.  
    Peter Boettcher, Aug 29, 2005
    #4
  5.  
    Roger Stafford, Aug 29, 2005
    #5
  6. <22-bradley.dialup.earth
    link.net>,
    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
    #6
  7. Christopher Hulbert

    Sui Huang Guest

    Hi! Roger:
    Thank you for your hard testing. Sorry about the late reply.
    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.)

    Thank you for your suggestion!

    Sui
     
    Sui Huang, Sep 1, 2005
    #7
  8. Christopher Hulbert

    Sui Huang Guest

    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
    #8
  9. Christopher Hulbert

    Bobby Cheng Guest

    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
    #9
  10. 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.
     
    Peter Boettcher, Sep 1, 2005
    #10
  11. Christopher Hulbert

    Sui Huang Guest

    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
    #11
  12. Christopher Hulbert

    Bobby Cheng Guest

    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.

    ---Bob.
     
    Bobby Cheng, Sep 1, 2005
    #12
  13. [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
    #13
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.