# How to obtain a Symplectic 4×4 matrix?

I have a problem in obtaining a $2n \times 2n$ Symplectic matrix $T$, with $n=2$. I couldn't find a direct command in Mathematica to achieve this.

Conditions:

Transpose[T].HH.T={{v1,0,0,0},{0,v2,0,0},{0,0,v1,0},{0,0,0,v2}
Transpose[T].JJ.T = JJ


where:

JJ={{0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 0}};
T={{T11,T12,T13,T14},{T21,T22,T23,T24},{T31,T32,T33,T34},{T41,T42,T43,T44}};
HH={{HH11,HH12,HH13,HH14},{HH21,HH22,HH23,HH24},{HH31,HH32,HH33,HH34},{HH41,HH42,HH43,HH44}};

• Have you looked at this article? library.wolfram.com/infocenter/MathSource/4779 – bill s Jul 14 '13 at 15:30
• I have already downloaded the code, but I don't know how to handle with it. Code is generalized and very complex so I can not use it in obtaining symplectic matrix. – Pipe Jul 14 '13 at 17:03
• @bill s can you halp me how to use this pakage to obtain symplectic matrix with condition. – Pipe Jul 17 '13 at 11:45
• that's a pretty general question. How about if you look at the code and try and figure it out, then pose a question on this site when you get stuck? – bill s Jul 17 '13 at 12:28
• thank you Bill, the problem is after running the package there many messages errors, mistakes in original code – Pipe Jul 17 '13 at 14:59

$HH$ appears linearly in $T^{\mathsf{T}}.HH.T=V$, and can be computed for a given symplectic matrix $T$ as $H=\left(T^{\mathsf{T}}\right)^{-1}.V.T^{-1}$.

solsH = Inverse[
Transpose[T]].{{v1, 0, 0, 0}, {0, v2, 0, 0}, {0, 0, v1, 0},
{0, 0, 0, v2}}.Inverse[T] // Simplify;


The symplectic matrix is not unique, so I am just going to get one that does not blow up $HH$. It turns out that in the expression for $solsH$ all the denominators are the same, so I don't have to manipulate that expression any more.

Flatten@Map[Denominator, solsH, {2}];
Differences[%]
(* {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} *)


The expression for the denominator that should not be zero.

den = Denominator[solsH[[1, 1]]];


Now I find an instance of $T$ such that $T^{\mathsf{T}}.JJ.T=JJ$ and $den\neq 0$.

solsT = FindInstance[
Join[Thread[Flatten[Transpose[T].JJ.T - JJ] == 0], {den != 0}],
Flatten[T]];


A set of solutions for $T$ and $HH$.

T /. solsT[[1]] // MatrixForm


$\left( \begin{array}{cccc} 1 & -2 & \frac{11}{2} & 13 \\ -1 & -\frac{1}{2} & 6 & 9 \\ 2 & -2 & 2 & 8 \\ -2 & 1 & 2 & 0 \\ \end{array} \right)$

solsH /. solsT[[1]] // MatrixForm


$\left( \begin{array}{cccc} 8 \text{v1}+68 \text{v2} & -2 \text{v2} & -13 \text{v1}-108 \text{v2} & -10 \text{v1}-73 \text{v2} \\ -2 \text{v2} & 8 \text{v1}+\text{v2} & 2 \text{v2}-9 \text{v1} & \frac{\text{v2}}{2}-14 \text{v1} \\ -13 \text{v1}-108 \text{v2} & 2 \text{v2}-9 \text{v1} & \frac{125 \text{v1}}{4}+173 \text{v2} & 32 \text{v1}+118 \text{v2} \\ -10 \text{v1}-73 \text{v2} & \frac{\text{v2}}{2}-14 \text{v1} & 32 \text{v1}+118 \text{v2} & 37 \text{v1}+\frac{325 \text{v2}}{4} \\ \end{array} \right)$

The two conditions are indeed satisfied.

Transpose[T].solsH.T /. solsT[[1]] // Simplify
(* {{v1, 0, 0, 0}, {0, v2, 0, 0}, {0, 0, v1, 0}, {0, 0, 0, v2}} *)

Transpose[T].JJ.T /. solsT[[1]]
(* {{0, 0, 1, 0}, {0, 0, 0, 1}, {-1, 0, 0, 0}, {0, -1, 0, 0}} *)

• So helpful, full answer with additional explanation of code. Thank you very much Suba, congrats. – Pipe Jul 14 '13 at 23:11
• I didn't mention that matrix HH already exist. So I need to obtain unique matrix T in function of matrix HH? To change solsT? T = Inverse[Transpose[T].H].V – Pipe Jul 15 '13 at 0:22
• @Pipe, if $HH$ is known, then $T^{\mathsf{T}}.\text{HH}.T=V$ gives 16 equations and there are 16 variables in $T$. So we will not be able to impose the condition that $T$ be symplectic. – Suba Thomas Jul 15 '13 at 3:59
• just a second, T is not unique, so if it is not unique, can I find one in agree with HH. Please, take a look in the beginning of the post, for special numerical case of given HH, can I obtain T? – Pipe Jul 15 '13 at 11:54
• The symplectic matrix is not unique. That is $T^{\mathsf{T}}.JJ.T=JJ$ has many solutions for $T$. But if we start from $T^{\mathsf{T}}.HH.T=V$ it appears that it will dictate the solution for $T$ and there seems to be no way to specify that it be symplectic. I will take another shot at your modified question later, assuming that you or someone else has not already figured it out. – Suba Thomas Jul 15 '13 at 16:15

I'll post an answer to an extremely closely related problem which I took some time sorting out today.

If I have some $$2n \times 2n$$ matrix $$J = \begin{pmatrix} 0 & - I \\ I & 0 \end{pmatrix}$$

    J = KroneckerProduct[{{0,-1},{1,0}},IdentityMatrix[n]]


and a second $$2n \times 2n$$ matrix $$H$$ which satisfies $$J H J = -H^*$$ for example

    H = RandomVariate[GaussianSymplecticMatrixDistribution[n]]


we know that $$H$$ is diagonalised by a symplectic unitary, that is $$H = U^\dagger S U \qquad \text{where} \qquad U^T J U = J.$$ and where $$S$$ is a diagonal matrix. However the unitaries spat out by eigensystem do not have this property. That is for

{S1, V} = Conjugate@Eigensystem[H];


we indeed obtain $$H = V^\dagger S V$$ but $$V$$ is not symplectic. This is because $$H$$ has a set of doubley degenerate subspaces and Eigensystem spits out some arbitrary basis vectors for these subpaces. This is remedied by using

    {S1, V} = ConjugateTranspose@SortBy[Transpose[Eigensystem[A]], Abs];

U = Join[#[[1 ;; -1 ;; 2]], #[[2 ;; -1 ;; 2]]] &@(Riffle[Table[1, d],
Conjugate@Table[V[[2 i]].J.V[[2 i - 1]], {i, 1, d}]] V);

S2 = Join[#[[1 ;; -1 ;; 2]], #[[2 ;; -1 ;; 2]]] &@S1;



which it can be verified satisfy the desired relations by checking

    (( ConjugateTranspose[U].DiagonalMatrix[S2].U - H ) // Norm // Chop ) == 0

(( Transpose[U].J.U - J ) // Norm // Chop ) == 0


return true.

• "However the unitaries spat out by Eigensystem do not have this property." - indeed, this is why there are specialized algorithms like the $SR$ algorithm for computing eigendecompositions; normal methods usually do not respect the symplectic structure. – J. M.'s technical difficulties May 19 at 2:33
• That's cool. Do you know if there's a option or something to get Eigensystem to do that? I couldn't find one, but in the end it seems little computational overhead is required to recover the symplectic structure from what mathematica does spit out. – ComptonScattering May 19 at 2:38
• No, unfortunately Eigensystem[] only does QR. So it comes down to a choice between post-processing the results of Eigensystem[] (which admittedly can still lead to some inaccuracy), or writing an implementation of $SR$ from scratch. (I'm not current on the literature anymore, so it is possible the method in that older article has been replaced by newer, more efficient ones. But it still does better than QR.) – J. M.'s technical difficulties May 19 at 2:42
• Useful for those willing to invest time finding a more precise solution than my quick fix! – ComptonScattering May 19 at 16:24