Fundamental matrix from camera matrices satisfying transpose property

Question

Assume we are given a pair of calibrated cameras with projection matrices $\mathbf \Pi_0$ and $\mathbf \Pi_1$. From Hartley and Zisserman (pg 246), the fundamental matrix is \begin{align} \mathbf F &= [\mathbf e_1]_\times \mathbf\Pi_1 \mathbf\Pi_0^\dagger \\ &= [\mathbf \Pi_1 \mathbf c_0]_\times \mathbf\Pi_1 \mathbf\Pi_0^\dagger \\ &= [-\mathbf \Pi_1 \mathbf M^{-1} \mathbf p_4]_\times \mathbf\Pi_1 \mathbf\Pi_0^\dagger \,, \end{align} where $\mathbf\Pi_0 = [\mathbf M \,|\, \mathbf p_4]$.

If I swap $\mathbf\Pi_0$ and $\mathbf\Pi_1$ (i.e., $\mathbf F' = [\mathbf e_0]_\times \mathbf\Pi_0 \mathbf\Pi_1^\dagger$), $\mathbf F'$ should satisfy the transpose property $\mathbf F' = \mathbf F^T$. I observe this property holds when I compute $\mathbf F, \mathbf F'$ using the 8-point algorithm, however, when I use the formula from Harley and Zisserman, the transpose property doesn't hold. Does anyone know what I'm doing wrong / if I need to normalize the projection matrices in some way?

Here is my implementation

import torch


def skew(x):
    return torch.tensor(
        [
            [0, -x[2], x[1]],
            [x[2], 0, -x[0]],
            [-x[1], x[0], 0],
        ]
    )


def compute_camera_center(P):
    """Hartley and Zisserman, pg 158"""
    return -P[:3, :3].inverse() @ P[:, 3]


def F(P0, P1):
    """Hartley and Zisserman, pg 246"""
    c0 = compute_camera_center(P0)
    e1 = P1 @ torch.concat([c0, torch.ones(1)])
    return skew(e1) @ P1 @ P0.pinverse()

For the example projection matrices below, the code produces valid fundamental matrices (i.e., they obey $\mathbf x_1^T \mathbf F \mathbf x_0 = 0$, but they are not transposes of each other

P0 = torch.tensor(
    [
        [-7.05e03, 3.45e03, 9.85e00, -4.90e05],
        [3.11e02, 9.13e02, 7.77e03, -1.57e06],
        [3.12e-01, 9.50e-01, -1.29e-03, -8.96e02],
    ]
)
P1 = torch.tensor(
    [
        [-1.33e03, -8.09e03, 2.89e01, 1.54e05],
        [-9.49e02, 4.71e01, 8.11e03, -1.38e06],
        [-1.00e00, 1.93e-02, 6.25e-05, -6.66e02],
    ]
)

print(F(P0, P1))
# tensor([[ 4.0068e+00,  1.0309e+03, -9.7512e+05],
#         [ 1.0242e+03, -4.1002e+00,  4.8394e+06],
#         [-9.8821e+05, -7.3285e+06,  2.4011e+09]])

print(F(P1, P0))
# tensor([[ 3.6779e+00,  9.4012e+02, -9.0703e+05],
#         [ 9.4625e+02, -3.7765e+00, -6.7266e+06],
#         [-8.9503e+05,  4.4420e+06,  2.2039e+09]])