Can a certain pair of expressions be compressed into one?

Question

I've been employing this pair of statements

z = Partition[Y1.Transpose[Y1], {2, 2}]; 
If[
  PositiveDefiniteMatrixQ[
    ArrayFlatten @ {{z[[1, 1]], z[[2, 1]]}, {z[[1, 2]], z[[2, 2]]}}] 
  == True, 
  r = r + 1]; 

in a larger program, where Y1 is a $4 \times 4$ matrix.

I'm testing the positive-definiteness of the "partial transpose" of z, in line with the quantum-information-theoretic Peres-Horodecki test for separability.

Can I eliminate the use of the variable z?
Can this task be accomplished in one command?
Can Flatten be employed for such a purpose?

Presumably, such a compression would lead to somewhat faster execution.

In light of the very detailed answer of Henrik Schumacher, I'm now including the larger program in which the pair of statements was embedded. He remarked: "Notice that a main part of the speed up comes from processing many matrices at once." I'd don't know whether such "processing" is possible/feasible in the framework of the larger program, since the generation of the matrices Y1 involves a number of other matrices, as well as a Check for an Orthogonalize error (for which, if it occurs I redo the calculation with added precision).

sp2 = x /. Solve[x^(37) == x + 1, x][[1]]; G = Array[1, 36]; Do[
 G[[i]] = N[i95/sp2^i], {i, 1, 36}]; rB = 0; hw = TimeUsed[]; Do[
 P = InverseCDF[NormalDistribution[0, 1], FractionalPart[G]]; 
 Label[wuh]; 
 Y1 = Check[(Orthogonalize[ArrayReshape[Take[P, {1, 16}], {4, 4}]] + 
      IdentityMatrix[4]).(ArrayReshape[Take[P, {17, 36}], {4, 5}]), 
   err; G1 = Array[1, 36]; Do[G1[[i]] = N[i95/sp2^i, 32], {i, 1, 36}];
    P = InverseCDF[NormalDistribution[0, 1], FractionalPart[G1]]; 
   Goto[wuh]]; z = Partition[Y1.Transpose[Y1], {2, 2}]; 
 If[PositiveDefiniteMatrixQ[
    ArrayFlatten@{{z[[1, 1]], z[[2, 1]]}, {z[[1, 2]], z[[2, 2]]}}] == 
   True, rB = rB + 1]; 
 If[Mod[i95, 1000000] == 0, 
  Print[{TimeUsed[] - hw, i95, rB, rB/i95, N[rB/i95, 20]}]; 
  hw = TimeUsed[]], {i95, 1, 100}]

Answer 1

Can Flatten be employed for such a purpose?

Yes. In this case, Flatten[#, {{1, 3}, {2, 4}}] & is equivalent to ArrayFlatten, but it really does not make a difference in runtime.

Can I eliminate the use of the variable z?

A somewhat more compact form can be obtained with the following function:

f[Y_] := ArrayFlatten[Transpose[Partition[Y.Transpose[Y], {2, 2}]]];

You can speed it up considerably by compiling it:

Block[{Y, YY},
  YY = Array[Compile`GetElement[Y, ##] &, {4, 4}];
  cf = With[{code = f[YY]},
    Compile[{{Y, _Real, 2}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
  ];

A timing test:

YY = RandomReal[{-1, 1}, {100000, 4, 4}];

a = f /@ YY; // RepeatedTiming // First
b = cf[YY]; // RepeatedTiming // First
a == b

0.881

0.0063

True

Notice that a main part of the speed up comes from processing many matrices at once.

This is how a refactorization of your code could look like:

r1 = Total[Boole[PositiveDefiniteMatrixQ /@ cf[YY]]]; // 
  AbsoluteTiming // First

0.232068

Ans a comparison to the original one:

r = 0;
Do[
    z = Partition[Y1.Transpose[Y1], {2, 2}];
    If[
     PositiveDefiniteMatrixQ[
      ArrayFlatten@{{z[[1, 1]], z[[2, 1]]}, {z[[1, 2]], z[[2, 2]]}}]
     ,
     r = r + 1
     ];
    , {Y1, YY}]; // AbsoluteTiming // First

r == r1

1.38389

True

Compiling also the computation of the eigenvalues ($4 \times 4$ matrices are the largest ones for which you can find closed form expressions for all eigenvalues), this can be further accelerated:

Block[{Y, YY},
  YY = Array[Compile`GetElement[Y, ##] &, {4, 4}];
  eigs4 = With[{code =
      N[Eigenvalues[YY, "Quartics" -> True]]
     },
    Compile[{{Y, _Complex, 2}},
     code,
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
  ];

Test:

ϵ = 1. 10^-12;
r2 = Total[
     UnitStep[
      UnitStep[Re[eigs4[cf[YY]] - ϵ]].ConstantArray[1, 4] - 
       4]]; // RepeatedTiming // First
r == r2

0.00037

True