Odd behaviour of 2x2-matrix and 2x1-vector operation

Consider a 2x2 matrix A with columns uncorrelated, and a vector x:

   A = [ a  b ]
       [-b' a']
   x = [ a']
       [ b']

where a and b are complex scalar. The product of A*x is supposed to be mathematically:

   A*x = [  a*a'+b*b' ] = [ a*a'+b*b' ]
         [-b'*a'+a'*b']   [ 0         ]

In matlab 2012b, I tried to compare A*x matrix operation and the element-wise operations such as a*a'+b*b' and -b'*a'+a'*b'. Both results happen to be different with some random pairs of a and b. a*a'*b*b' always shows positive real, but matrix operation shows some residual on the imaginary part sometime. Even -b'*a'+a'*b' entry of matrix operation shows such erroneous imaginary part, which is supposed to be zero!

The following script finds the mismatch result and shows the hex forms of the results.

y0, y1, y2, z show the result of different operations:

• y0 = A*x; % matlab matrix multiply script
• y1 = mtimes(A,x); % mtimes operation
• y2 = mvmult_lapack (A, x); % private code using lapack routine (ZGEMV)
• z = [a*a'+b*b';-b'*a'+a'*b']; % element-wise operation

Result : y2 and z are always same and real, but y0 and y1 have erroneous imaginary part. Note that the real parts of y0 and y1 are different from the real parts of y2 and z by eps, say a difference of the last bit of floating-point format.

Discuss and Question : It is not clear why the A*x script shows the difference from the element-wise operation. It could not be a problem with LAPACK/BLAS. There seems parser or storage register problem in machine code when the matrix operation is called. Is this only preblem with 2012b or my machine?

Test script :

randn('seed',0)
ii=0;
while 1,
    ii = ii + 1;
    a = randn+i*randn;
    b = randn+i*randn;
    %
    A = [a b; -b' a'];
    x = [a'; b'];
    y0 = A*x;
    y1 = mtimes(A,x);
    y2 = mvmult_lapack (A, x);
    z = [a*a'+b*b';-b'*a'+a'*b'];
    if sum(y0~=z) % break when the results are different!
        break
    end
    if ii==10000
        break
    end
end
ii
%
a
b
y0
y1
y2
z
y0==z
[num2hex(real(y0)), repmat(' ',[2 1]), num2hex(imag(y0))]
[num2hex(real(y1)), repmat(' ',[2 1]), num2hex(imag(y1))]
[num2hex(real(y2)), repmat(' ',[2 1]), num2hex(imag(y2))]
[num2hex(real(z)), repmat(' ',[2 1]), num2hex(imag(z))]

Output :

 ii =

      1

 a =

            1.16495351050066 +     0.626839082632431i 

 b =

          0.0750801546776829 +     0.351606902768522i

 y0 = 

            1.87930836084417 -  3.46944695195361e-18i
                           0 -  2.08166817117217e-17i

 y1 =

            1.87930836084417 -  3.46944695195361e-18i
                           0 -  2.08166817117217e-17i

 y2 =

           1.87930836084417
                          0

 z =

           1.87930836084417
                          0

 ans =

      0
      0

 ans =

 3ffe11a5a4cecd1e bc50000000000000
 0000000000000000 bc78000000000000

 ans =

 3ffe11a5a4cecd1e bc50000000000000
 0000000000000000 bc78000000000000

 ans =

 3ffe11a5a4cecd1d 0000000000000000
 0000000000000000 0000000000000000

 ans =

 3ffe11a5a4cecd1d 0000000000000000
 0000000000000000 0000000000000000

[Stephen23]

Maybe you should start a Newsgroup thread about this... it seems like an interesting issue to discuss.

Answer 1

[Roger Stafford]

Both the differences in real parts and the spurious imaginary parts in y0 and y1 are extremely small and are undoubtedly due to rounding differences in the four different methods of computation. Your results show that the real parts differ only in the least significant bits and the imaginary parts are likewise extremely small.

You should realize that simply rearranging the ordering of computations can result in slightly differing results. An example I like to give is:

   format hex
   (3/14+3/14)+15/14 --> 3ff8000000000000
   3/14+(3/14+15/14) --> 3ff7ffffffffffff

Although both should have 1.5 as a result, the second one comes out one bit different from that. That is because there are two operations performed with a rounding after each one, and this is done in a different order in the two cases with a resulting difference in the final result.

The same applies to the operations performed in your four methods. The sequencing of operations is different in the four methods and the results therefore don't exactly match. This is something that has to be accepted in performing numerical computations using digital computers with a finite number of bits accuracy.

[Sung-Eun Jo]

Speaking of computation ordering and rounding-off, I thought that 2x2 complex matrix operation is concerning only 2 terms each entry such as a*a'+b*b' and -b'*a'+a'*b'. I was wonder how the ordering problem takes place in such 2 term case.

After your comment, I realized that 2 term complex operation is extended to 4 term real operation each real and imaginary part:

 a*a'+b*b'
 = (ar+i*ai)*(ar-i*ai)+(br+i*bi)*(br-i*bi)
 = (ar*ar + ai*ai + br*br + bi*bi) + i*(ai*ar - ar*ai + bi*br - br*bi)

-b'*a'+a'*b'
= (-br+i*bi)*(ar-i*ai)+(ar-i*ai)*(br-i*bi)
= (-br*ar + bi*ai + ar*br - ai*bi) + i*(br*ai + bi*ar - ar*bi - ai*br)

 ar = real(a)
 ai = imag(a)
 br = real(b)
 bi = imag(b)

 zr = -br*ar + bi*ai + ar*br - ai*bi
 zi = br*ai + bi*ar - ar*bi - ai*br

 zr2 = (-br*ar + bi*ai) + (ar*br - ai*bi)
 zi2 = (br*ai + bi*ar) + (- ar*bi - ai*br)

 zi =
     -2.08166817117217e-17
 zi2 =
     0

In the above expansion, I noticed two feasible ordering as in zi and zi2. After comparing the results of zi and zi2, we prefer the computation order same as in zi2, which is meaningful when we are considering matrix entry-wise operation.

However, it is still hard to explain WHY the computation ordering is different between the direct computations [a*a'+b*b'; -b'*a'+a'*b'] and the matrix-vector operation A*x.

The result shows that the direct computation order (and Lapack zgemv as well) provides better round-off performance. If so, we may conclude that the complex matrix product A*x in matlab seems to do some weird computation ordering.

Lapack zgemv routine seems to use the same order of the direct computation: It is because the matrix complex entry-product takes place first and entry-sum later:

(-b'*a') + (a'*b')
= ( (-br*ar + bi*ai) + (ar*br - ai*bi) ) + i*( (br*ai + bi*ar) + (- ar*bi - ai*br) )

I guess that the complex matrix-vector product, A*x, should use the same computation order as in Lapack.

-- jo

[Roger Stafford]

The ordering (grouping) for the imaginary part

   ( (br*ai) + (bi*ar) ) + ( (-ar*bi) + (-ai*br) )

should always produce an exact zero in matlab. That is because two-operand addition and multiplication are always commutative in the IEEE 754 standard. However the ordering

   ( ( (br*ai) + (bi*ar) ) + (-ar*bi) ) + (-ai*br)

or other similar orderings could easily yield a nonzero, (though of course an extremely small one,) because the associative law does not always hold. Something like the latter method is therefore probably what 'mtimes' is using.

[Sung-Eun Jo]

I agree. In matlab, 'mtimes' and * operator for complex data should be careful. Their computation ordering may cause numerical errors even in trivial cases such as perfect zero or real value expected.