Givens rotation: What is the Mathematica equivalent to MATLAB's planerot function?

Question

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.

Answer 1

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.

Answer 2

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