2x2 orthogonal matrix multiplication error

[s.jo]

Hello all,

In octave(3.8.2), 2x2 orthogonal matrix M shows numerical error after
multiplied by own transpose.

Say, 
   M*M' is not equal to eye(2)!

There are errors in off-diagonal entries.
But, if we try it with element-wise multiplication, it works well.

Say, 
Mt = M';
[Mt(1, :) * M(:, 1), Mt(1, :) * M(:, 2); Mt(2, :) * M(:, 1), Mt(2, :) * M(:,
2)] is equal to eye(2)!!

I check this also with matlab, but it works well.
I suspect that octave has some unreliable routine with
matrix-multiplication.
Or, is this only problem with my cygwin machine?

I also try it with the complex rotation matrix (say Hermittian 2x2 matrix).
It also fails.
I generate rotation matrix with cos and sin functions:

M = [cs sn; -sn' cs'] 

The following script shows these result.

Such matrix orthogonality easily fails with large size matrix due to
numerical errors.
But 2x2 matrix is so small. It should work always with such orthogonality
like matlab do, I guess.

Can somebody help this how it figure out?

Many thanks,

jo

octave> uname
ans =

  scalar structure containing the fields:

    sysname = CYGWIN_NT-6.1
    nodename = mkkim
    release = 1.7.32(0.274/5/3)
    version = 2014-08-13 23:03
    machine = i686

octave> test_complex_mat_mult
+ ## Test of rotation matrix multiplication
+ ## Error in matrix multiplication
+ ## There is some error with column-uncorrelated matrix multiplication.
+ version
ans = 3.8.2
+ ## Complex rotation matrix
+ ## Generate a complex rotation matrix
+ theta = randn * pi;
+ cs = cos (theta)
cs =   -0.56511
+ phase = sign (randn + i * randn);
+ sn = sin (theta) * phase
sn =   -0.39729 -    0.72306i
+ M = [cs, sn; -sn', cs]
M =

    -0.56511 +          0i    -0.39729 -    0.72306i
     0.39729 -    0.72306i    -0.56511 +          0i

+ ##
+ ## Check 
+ abs (cs) ^ 2 + abs (sn) ^ 2 == 1
ans =          0
+ ##
+ ## Matrix products: M'*M, M*M'
+ M' * M
ans =

           1 +          0i  -9.8256e-18 - 2.6834e-18i
  -9.8256e-18 + 2.6834e-18i           1 +          0i

+ M * M'
ans =

           1 +          0i  9.8256e-18 + 2.6834e-18i
  9.8256e-18 - 2.6834e-18i           1 +          0i

+ ##
+ ## Element-wise products
+ Mt = M';
+ M (1, :) * Mt
ans =

           1           0

+ M (2, :) * Mt
ans =

           0           1

+ Mt * M (:, 1)
ans =

           1 +          0i
  -9.8256e-18 + 2.6834e-18i

+ Mt * M (:, 2)
ans =

  -9.8256e-18 - 2.6834e-18i
           1 +          0i

+ [Mt(1, :) * M(:, 1), Mt(1, :) * M(:, 2); Mt(2, :) * M(:, 1), Mt(2, :) *
M(:, 2)]
ans =

           1           0
           0           1

+ ## Real rotaion matrix
+ ## Generate a real rotation matrix
+ theta = randn * pi;
+ cs = cos (theta)
cs =    0.55482
+ phase = 1;
+ sn = sin (theta) * phase
sn =   -0.83197
+ M = [cs, sn; -sn', cs]
M =

     0.55482    -0.83197
     0.83197     0.55482

+ ##
+ ## Check 
+ abs (cs) ^ 2 + abs (sn) ^ 2 == 1
ans =          1
+ ##
+ ## Matrix products: M'*M, M*M'
+ M' * M
ans =

           1  -1.0137e-17
  -1.0137e-17           1

+ M * M'
ans =

           1  1.0137e-17
  1.0137e-17           1

+ ##
+ ## Element-wise products
+ Mt = M';
+ M (1, :) * Mt
ans =

           1           0

+ M (2, :) * Mt
ans =

           0           1

+ Mt * M (:, 1)
ans =

           1
  -1.0137e-17

+ Mt * M (:, 2)
ans =

  -1.0137e-17
           1

+ [Mt(1, :) * M(:, 1), Mt(1, :) * M(:, 2); Mt(2, :) * M(:, 1), Mt(2, :) *
M(:, 2)]
ans =

           1           0
           0           1

+ ## Note
+ ## M' seems to have unwanted numerical error in octave!
+ ## The element-wise product is good, but matrix product has some error 
+ ## where the entries should be zeros.
+ ## If we try real rotation matrix, it shows same.
+ endscript;




--
View this message in context: 
http://octave.1599824.n4.nabble.com/2x2-orthogonal-matrix-multiplication-error-tp4667797.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.