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.
planerot
gives the same output as that of your Mathematica implementation: octave-online.net/… $\endgroup$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
andConjugate[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$