# Can one use the new random matrices features of Version 11 in addressing a certain octonionic-based question?

I have cross-posted the question

https://math.stackexchange.com/questions/2012301/can-one-test-an-octonionic-interpretation-for-a-conjecture-apparently-valid-in

and

https://physics.stackexchange.com/questions/292540/can-one-test-an-octonionic-interpretation-for-a-quantum-information-conjecture

It seems feasible that the new random matrices features in Version 11 might facilitate addressing this question. Any specific thoughts in this matter? One of the new commands is CircularQuaternionMatrixDistribution. Can one create an octonionic counterpart?

The two formulas $P_1(\alpha)$ and $P_2(\alpha)$ were developed based solely on analyses of matrices with real and complex (and not quaternionic and octonionic) entries. To be more specific, the ascending moments of determinants (emphasis added) of the 4 x 4 matrices and of their “partial transposes” were computed, and formulas found for them. (These were, then, used in the Mathematica density approximation procedure of Provost [http://www.mathematica-journal.com/issue/v9i4/contents/DensityApproximants/DensityApproximants.pdf], to eventually arrive at the expressions for $P_1(\alpha)$ and $P_2(\alpha)$ .)

The two formulas (Charles Dunkl observed) could be absorbed into one, by regarding the parameter in the complex case to be twice that in the real case (hence the apparent [Dyson-index-like] connection to random matrix theory).

Now, although the calculation of determinants is straightforward with matrices the entries of which are restricted to real and complex values, it becomes more subtle with the quaternions, and, a fortiori, it would seem with the octonions. E. H. Moore (Bull. Amer. Math. Soc. 28 [1922], 161-162) gave a definition in the quaternionic case—and Wikipedia has a brief article, I see, about the “Dieudonne determinant” (“which is a generalization of the determinant of a matrix over division rings and local rings”). So, I think the originally stated problem I posed hinges on to what extent the formula Dunkl developed can be “extrapolated” to the octonionic domain. (I note, however, that Fei and Joynt in the cited paper appear to have by-passed the use of determinants, in their quaternionic analysis).

In preparing this “answer”, I found a (rather remarkable) series of (unpublished) June 2012 emails from Dunkl in which he does a highly in-depth [using Maple] analysis of the use of the Moore determinant in the quaternionic case, apparently succeeding in confirming its appropriateness. Here is part of his treatment (this, of course, deal with the quaternionic scenario, and the original octonionic question remains). Dunkl writes more of interest in this series of detailed emails (but I don’t see how to really present his remarks here).

The Maple code Dunkl employed for the application of the Moore determinant to the quaternionic case was:

qm := proc (z1, z2) local zq1, zq2, w1, w2, w3, w4; global qco, iq, jq, kq; zq1 := qco(z1); zq2 := qco(z2); w1 := zq1[1]*zq2[1]-zq1[2]*zq2[2]-zq1[3]*zq2[3 ]-zq1[4]*zq2[4]; w2 := zq1[1]*zq2[2]+zq1[2]*zq2[1]+zq1[3]*zq2[4]-zq1[4]*zq2[3] ; w3 := zq1[1]*zq2[3]+zq1[3]*zq2[1]+zq1[4]*zq2[2]-zq1[2]zq2[4]; w4 := zq1[1] zq2[4]+zq1[4]*zq2[1]+zq1[2]*zq2[3]-zq1[3]*zq2[2]; w1+w2*iq+w3*jq+w4*kq end proc; qconj := proc (f) options operator, arrow; subs({kq = -kq, iq = -iq, jq = -jq} ,f) end proc; qco := proc (f) local fq; global iq, jq, kq, lq; fq := collect(f,[iq, jq, kq]) ; lq[2] := coeff(fq,iq); lq[3] := coeff(fq,jq); lq[4] := coeff(fq,kq); lq[1] := subs(iq = 0,jq = 0,kq = 0,f); [lq[1], lq[2], lq[3], lq[4]] end proc; qdet4x := [[[1, 1], [2, 2], [3, 3], [4, 4], 1], [[4, 4], [3, 3], [1, 2], [2, 1 ], -1], [[4, 4], [2, 3], [3, 2], [1, 1], -1], [[4, 4], [2, 2], [1, 3], [3, 1], -1], [[3, 3], [2, 4], [4, 2], [1, 1], -1], [[3, 4], [4, 3], [2, 2], [1, 1], -1 ], [[3, 4], [4, 3], [1, 2], [2, 1], 1], [[2, 3], [3, 2], [1, 4], [4, 1], 1], [ [2, 4], [4, 2], [1, 3], [3, 1], 1], [[4, 4], [1, 2], [2, 3], [3, 1], -1], [[4, 4], [1, 3], [3, 2], [2, 1], -1], [[3, 3], [1, 4], [4, 2], [2, 1], -1], [[2, 2] , [1, 3], [3, 4], [4, 1], -1], [[2, 2], [1, 4], [4, 3], [3, 1], -1], [[2, 3], [3, 4], [4, 2], [1, 1], -1], [[2, 4], [4, 3], [3, 2], [1, 1], -1], [[1, 2], [2 , 3], [3, 4], [4, 1], 1], [[1, 2], [2, 4], [4, 3], [3, 1], 1], [[1, 3], [3, 2] , [2, 4], [4, 1], 1], [[1, 3], [3, 4], [4, 2], [2, 1], 1], [[1, 4], [4, 2], [2 , 3], [3, 1], 1], [[1, 4], [4, 3], [3, 2], [2, 1], 1], [[3, 3], [2, 2], [1, 4] , [4, 1], -1], [[3, 3], [1, 2], [2, 4], [4, 1], -1]]; qmdet4 := proc (mx) local dt, i, tm, ppq; global qm, qdet4x; dt := 0; for i to 24 do tm := op(i,qdet4x); ppq := qm(mx[tm[1][1],tm[1][2]],qm(mx[tm[2][1],tm[2] [2]],qm(mx[tm[3][1],tm[3][2]],mx[tm[4][1],tm[4][2]]))); dt := dt+tm[5]* simplify(ppq) end do; simplify(dt) end proc;

Also, here is the list of the 24 (4!) factors of the Moore determinant, in order, with the sign:

# this is a list of the 24 factors, in order, with the sign

qdet4x;

[[[1, 1], [2, 2], [3, 3], [4, 4], 1],

    [[4, 4], [3, 3], [1, 2], [2, 1], -1],

[[4, 4], [2, 3], [3, 2], [1, 1], -1],

[[4, 4], [2, 2], [1, 3], [3, 1], -1],

[[3, 3], [2, 4], [4, 2], [1, 1], -1],

[[3, 4], [4, 3], [2, 2], [1, 1], -1],

[[3, 4], [4, 3], [1, 2], [2, 1], 1],

[[2, 3], [3, 2], [1, 4], [4, 1], 1],

[[2, 4], [4, 2], [1, 3], [3, 1], 1],

[[4, 4], [1, 2], [2, 3], [3, 1], -1],

[[4, 4], [1, 3], [3, 2], [2, 1], -1],

[[3, 3], [1, 4], [4, 2], [2, 1], -1],

[[2, 2], [1, 3], [3, 4], [4, 1], -1],

[[2, 2], [1, 4], [4, 3], [3, 1], -1],

[[2, 3], [3, 4], [4, 2], [1, 1], -1],

[[2, 4], [4, 3], [3, 2], [1, 1], -1],

[[1, 2], [2, 3], [3, 4], [4, 1], 1],

[[1, 2], [2, 4], [4, 3], [3, 1], 1],

[[1, 3], [3, 2], [2, 4], [4, 1], 1],

[[1, 3], [3, 4], [4, 2], [2, 1], 1],

[[1, 4], [4, 2], [2, 3], [3, 1], 1],

[[1, 4], [4, 3], [3, 2], [2, 1], 1],

[[3, 3], [2, 2], [1, 4], [4, 1], -1],

[[3, 3], [1, 2], [2, 4], [4, 1], -1]]


# this adds the 24 terms to get the Moore determinant

print(qmdet4);

proc(mx) local dt, i, tm, ppq; global qm, qdet4x; dt := 0; for i to 24 do tm := op(i, qdet4x); ppq := qm(mx[tm[1][1], tm[1][2]], qm( mx[tm[2][1], tm[2][2]], qm( mx[tm[3][1], tm[3][2]], mx[tm[4][1], tm[4][2]]))) ; dt := dt + tm[5]*simplify(ppq) end do; simplify(dt) end proc

# a14,a23,a41,a32

mtqp;

    [    2
[h[1]  , 0 , 0 ,

]
h[2] g2[0] + h[2] g2[1] iq + h[2] g2[2] jq + h[2] g2[3] kq]

[        2
[0 , h[2]  ,

h[1] g1[0] + h[1] g1[1] iq + h[1] g1[2] jq + h[1] g1[3] kq ,

]
0]

[
[0 ,

h[1] g1[0] - h[1] g1[1] iq - h[1] g1[2] jq - h[1] g1[3] kq ,

2        2        2        2       2    ]
g2[0]  + g2[1]  + g2[2]  + g2[3]  + h[3]  , 0]

[
[h[2] g2[0] - h[2] g2[1] iq - h[2] g2[2] jq - h[2] g2[3] kq ,

2        2        2        2       2]
0 , 0 , g1[0]  + g1[1]  + g1[2]  + g1[3]  + h[4] ]


# this is the PT det (Moore formula)

det4;

     2      2       2      2       2      2       2      2


-(-h[2] g2[0] - h[2] g2[1] - h[2] g2[3] - h[2] g2[2]

           2      2       2      2       2      2       2      2
+ h[1]  g1[0]  + h[1]  g1[1]  + h[1]  g1[2]  + h[1]  g1[3]

2     2       2      2       2      2       2      2
+ h[4]  h[1] ) (h[1]  g1[3]  + h[1]  g1[0]  + h[1]  g1[1]

2      2       2     2       2      2       2      2
+ h[1]  g1[2]  - h[2]  h[3]  - h[2]  g2[0]  - h[2]  g2[1]

2      2       2      2
- h[2]  g2[2]  - h[2]  g2[3] )


• Maybe you could describe which random matrix functions you mean, and how they might apply? – bill s Nov 14 '16 at 3:19
• Thanks, bill s. Peter Forrester posted (arxiv.org/abs/1610.08081) a paper, "Octonions in Random Matrix Theory", outlining certain difficulties in dealing with a theory of Hermitian random matrices with octonion entries. The case N=3 is problematical, while the conjecture in question deals with the even more challenging case N=4. However, there seems to must be some way of addressing the question that quite naturally arises from application of the formulas I give. So, I want/need 4 x 4 Hermitian random matrices with octonion entries (and unit trace). – Paul B. Slater Nov 14 '16 at 5:07
• What is the purpose of posting this and your self-answer on three different sites? Please format your code snippets correctly, and make sure that the code is correct and complete. – Yves Klett Nov 26 '16 at 22:16
• Thanks, Yves Klett. The problem seemed to have mathematical/physical/Mathematica aspects. Sorry, for possibly over-reaching here. – Paul B. Slater Nov 29 '16 at 12:33
• I have expanded my answer into "Octonionic separability probability conjectures" (arxiv.org/abs/1612.02798) – Paul B. Slater Dec 12 '16 at 3:19

It appears that we can make meaningful progress in addressing this problem--via numerical computations--but without yet fully resolving a number of issues.

The question as put pertains to $4 \times 4$ (density) matrices--and, in this regard, we seek to extend the quite recent analyses of $2 \times 2$ and $3 \times 3$ "Wishart matrices ($W$) with octonion entries" of Peter Forrester (sec. 3 of https://arxiv.org/pdf/1610.08081.pdf). He employed Cholesky decompositions $W=T^{\dagger} T$. Accordingly, we start with $4 \times 4$ null matrices $T$ and fill their six upper triangular off-diagonal entries with octonions, the eight independent components of each of the six distributed as standard Gaussians.

Next, we fill the four diagonal entries with values that are the square roots of Gamma distribution variates. For the $2 \times 2$ case, Forrester employs $\Gamma[a+1,2]$ and $\Gamma[a+5,2]$. In the $3 \times 3$ instance, he utilizes $\Gamma[a +4 (i-1),2]$, $i=1,2,3$. (The parameter $a$ appears to not need to be fully specified, but must be large enough that the Gamma distributions is well-defined.) For the $4 \times 4$ case, we initially employed $\Gamma[a+3 (i-1),2]$, $i=1,\ldots,4$ (but, admittedly, have no compelling argument for this particular choice).

To further proceed, we relied upon the suite of Mathematica programs made available by Tevian Dray and Corinne A. Manogue in their paper, "Finding octonionic eigenvectors using Mathematica" (Computer Physics Commun 115 [1998], 536-547). This allowed us (using, in particular, their commands "omult[p,q]" and "MMult[X,Y]") to generate random Wishart matrices $W=T^{\dagger} T$.

For each such matrix, we sought to compute its determinant and then test its positivity. For this purpose we attempted to employ "Theorem 5.3. (Laplace expansion)" in the 2010 paper of Jianquan Liao, Jinxun Wang and Xingmin Li in volume 20 of Anal. Theory Appl., entitled "The all-associativity of octonions and its applications"(http://link.springer.com/article/10.1007/s10496-010-0326-2). (We followed the "template" of the Laplace expansion of a $4 \times 4$ matrix by $2 \times 2$ "complementary minors" presented in https://en.wikipedia.org/wiki/Laplace_expansion, and utilized the Mathematica command "Odet[X]" in the Dray-Manogue package for the computation of the $2 \times 2$ minors.)

At this point, we were prepared for our simulation of the Wishart matrices $W=T^{\dagger} T$. At first, we set $a=1$. For one thousand such random matrices, we found (using the Laplace expansion algorithm) all of their determinants to be positive. (Let us note that Forrester in simulating ten thousand $3 \times 3$ Wishart matrices, found 5,500 of them to have negative determinants--but the particular value of $a$ employed was not indicated.)

Now, to address the underlying/central question of the value of the "[Hilbert-Schmidt] separability probability of two-qubit density matrices with octonionic entries", we (again, using the Laplace expansion routine) computed the determinants of the "partial transposes" of the one thousand (positive determinant) Wishart matrices. Of the one thousand partial transposes, 651 had positive determinants, giving us a separability probability of 0.651, orders of magnitude larger than the conjectured value (see initial question) of $P_1(4)=P_2(4) = \frac{44482}{4091349} \approx 0.0108722$.

But now, in our most interesting finding, we ascertained that one could readily "tune" the separability probability by the choice of the parameter $a$. So, for $a=\frac{1}{128}$, we obtained 996 Wishart matrices with positive determinants. Of these, only 13 now had "positive partial transposes (PPT's), giving us a separability probability of $\frac{13}{996} \approx 0.0130522$, much closer to the conjectured value.

So, at this point in time, the main question seems to be does there exist some underlying principle for the choice of the parameter $a$ (as well as what are the appropriate Gamma distributions to employ for setting the diagonal entries of the random upper triangular matrices $T$ in the Cholesky decompositions $W=T^{\dagger} T$ for the random $4 \times 4$ Wishart matrices)? Also, does it suffice to test positivity to only employ determinants, as opposed to the (more computationally problematical) eigenvalues?