3
$\begingroup$

Matlab has a function called planerot which takes a two-component column vector $x$ as input and returns a $2 \times 2$ matrix $G$ and a two component column vector $y$ such that $y = G.x$ and $y_2 = 0$. The documentation is available here.

I am trying to implement it using the ROTG routine of the BLAS package in Mathematica as follows

Needs["LinearAlgebra`BLAS`"];
planerot[Xcol_] := Module[{a, b, a0, b0, c, s, y, givens},
  {a, b} = {a0, b0} = Xcol;
  ROTG[a, b, c, s];
  givens = {{c, s}, {-Conjugate[s], c}};
  y = givens . Transpose[{a0, b0}];
  {givens, y}]

I have tested this using the example given here and it works BUT if I pass a complex vector as the argument, Mathematica and Matlab gives different answers.

First, using $x = (3 ,4)^T$, we can see that Mathematica gives,

In[ ]:= {G, y} = planerot[{3, 4}] // N

Out[ ]= {{{0.6, 0.8}, {-0.8, 0.6}}, {5., 0.}}

which matches with the values given by MATLAB.

But using an example from Mathematica's documentation, our tiny method will give the following output for $x = (-2+3i,-2i)^T$,

In[808]:= {Gt, y} = planerot[{-2. + 3. I, -2. I}] // Chop 

Out[808]= {{{0.874475, -0.403604 - 0.269069 I}, {0.403604 - 
    0.269069 I, 0.874475}}, {-2.28709 + 3.43063 I, 0}}

whereas, MATLAB gives us

>> [G,y]= planerot([-2+3i,-2i]')

G =

  -0.4851 + 0.7276i   0.0000 - 0.4851i
   0.0000 - 0.4851i  -0.4851 - 0.7276i


y =

    4.1231
         0

In both cases, $y_2 = 0$ but $y_1$ and the Givens matrix are different. How can I resolve this discrepancy? Or is there an open-source implementation of planerot that I can look at?

I would really appreciate your help.

$\endgroup$
3
  • 2
    $\begingroup$ Funny, at the moment Octave's planerot gives the same output as that of your Mathematica implementation: octave-online.net/… $\endgroup$
    – xzczd
    Feb 22, 2023 at 8:55
  • $\begingroup$ @xzczd I think the differences in the result arise from differences in definition. Although, thanks to the two answers here I can see how the differences arise, I am frankly a bit perplexed as to how to proceed. Which algorithm is "correct"? $\endgroup$
    – noir1993
    Feb 22, 2023 at 12:18
  • 2
    $\begingroup$ Hey, those solutions only differ by a multiplication with a complex rotation, i.e., with d = DiagonalMatrix[{Exp[I Arg[-2.28709 + 3.43063 I]], -I Exp[-I Arg[-0.403604-0.269069I]]}]. If you compute {{0.874475, -0.403604 - 0.269069 I}, {0.403604 - 0.269069 I, 0.874475}}.d and Conjugate[d].{-2.28709 + 3.43063 I, 0}, you'll get back MATLAB's solution. Now, which complex rotation is "correct" is a matter of taste. The same normalisation choice applies to eigenvalue decomposition, SVD, etc., where the solutions of Mathematica and MATLAB might vary. $\endgroup$ Feb 22, 2023 at 16:22

2 Answers 2

4
$\begingroup$

I don't know how this does make sense ($G$ is not orthogonal for the example complex input and $G.x \neq y$ !), but it seems to reproduce the results:

ClearAll[planerot];
planerot[x : {x1_, x2_}] :=
 With[{norm = Norm[x]},
  {{{x1, x2},
    Conjugate[{-x2, x1}]}/norm,
   {norm, 0}}]

Frankly the fact that this function reproduces two results in the question doesn't cast much confidence in me for this being actually correct solution.

Numerical stability is a question I won't touch with a ten-foot pole, maybe you should check out the Wikipedia entry for Givens rotation.

$\endgroup$
4
  • 2
    $\begingroup$ Your function produces a matrix $G$ whose $G_{44}$ element is off by a sign. MATLAB gives -0.4851 - 0.7276i and your MMA function gives -0.485071 + 0.727607 i. I suppose taking a complex conjugate will fix the problem. $\endgroup$
    – noir1993
    Feb 22, 2023 at 12:30
  • 1
    $\begingroup$ @noir1993 Indeed! I have no idea why conjugates make this "work", but after that fix the matrix is at least unitary. I'm just wondering how results in general compare with the Matlab implementation... $\endgroup$
    – kirma
    Feb 22, 2023 at 12:47
  • 2
    $\begingroup$ @noir1993 and kirma, I just noticed we've all missed a feature of MATLAB core language, see the updated answer of mine. $\endgroup$
    – xzczd
    Feb 22, 2023 at 13:02
  • $\begingroup$ @xzczd Thank you! Great catch. Now my code gives me the correct $y$ but the Givens matrix is the complex conjugate of the one given by MATLAB. $\endgroup$
    – noir1993
    Feb 22, 2023 at 13:25
4
$\begingroup$

Definitely not the most efficient, but the most straightforward solution in my mind:

planerot = x |-> With[{y = {Norm@x, 0}}, {RotationMatrix@{x, y}, y}];

planerot@{3, 4.}
(* {{{0.6, 0.8}, {-0.8, 0.6}}, {5., 0}} *)

planerot@Conjugate@{-2. + 3 I, -2 I} // Chop
(* {{{-0.485071 + 0.727607 I,  0.       - 0.485071 I}, 
     { 0.       - 0.485071 I, -0.485071 - 0.727607 I}}, 
    {4.12311, 0}} *)

Notice I've added Conjugate in 2nd example, because in MATLAB ' is the shorthand for ctranspose (complex conjugate transpose). BTW, shorthand for transpose is .'.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.