I cannot help but suspect that the working function you've provided is less efficient that it can be, but I must also admit that I have only dabbled in computer vision myself. Here's a suggestion for how it might be improved.

As you've mentioned, the fundamental matrix $F$ is defined as the matrix that satisfies $x_{1}^{T} F x_{2} = 0$ for all points $x_{1}$ and $x_{2}$, in the first and second images, respectively, which correspond to the same point $z$ on the 3D object of interest. That is, for all $z$, $F$ satisfies $$z^{T} P_{1}^{T} F P_{2} z = 0.$$

A quadratic form is constantly $0$ if and only if its matrix is skew-symmetric (See here for justification). Therefore $F$ is the matrix such that $$P_{1}^{T} F P_{2} + P_{2}^{T} F^{T} P_{1} = 0.$$ Applying the vectorization operator to both sides of this equation and using the compatibility of vectorization with kronecker product (see here for justification), we have $$(P_{2}^{T} \otimes P_{1}^{T})\, \text{vec}(F) + (P_{1}^{T} \otimes P^{T}_{2}) \, \text{vec}(F^{T}) = 0.$$ Denote by $K$ the commutation matrix of $F$'s dimension. Then $$(P_{2}^{T} \otimes P_{1}^{T})\, \text{vec}(F) + (P_{1}^{T} \otimes P^{T}_{2}) \, K^{T} \, \text{vec}(F) = 0.$$

Finally, $F$ can be solved for by unvectorizing the solution to $$\left((P_{2}^{T} \otimes P_{1}^{T}) + (P_{1}^{T} \otimes P^{T}_{2}) \, K^{T} \right) \, \text{vec}(F) = 0,$$ which you can solve with scipy. Of course, the commutation matrix should be computed once and used throughout the for loop. I prefer this approach because it works with the camera matrices directly. I suspect it will be faster but cannot confirm.

Note also that the fundamental matrix between cameras $i$ and $j$ is the transpose of the fundamental matrix between cameras $j$ and $i$. This allows you to cut the number of fundamental matrices you are computing in half.

I cannot help but suspect that the working function you've provided is less efficient than it can be, but I must also admit that I have only dabbled in computer vision myself. Here's a suggestion for how it might be improved.

As you've mentioned, the fundamental matrix $F$ is defined as the matrix that satisfies $x_{1}^{T} F x_{2} = 0$ for all points $x_{1}$ and $x_{2}$, in the first and second images, respectively, which correspond to the same point $z$ on the 3D object of interest. That is, for all $z$, $F$ satisfies $$z^{T} P_{1}^{T} F P_{2} z = 0.$$

A quadratic form is constantly $0$ if and only if its matrix is skew-symmetric (See here for justification). Therefore $F$ is the matrix such that $$P_{1}^{T} F P_{2} + P_{2}^{T} F^{T} P_{1} = 0.$$ Applying the vectorization operator to both sides of this equation and using the compatibility of vectorization with kronecker product (see here for justification), we have $$(P_{2}^{T} \otimes P_{1}^{T})\, \text{vec}(F) + (P_{1}^{T} \otimes P^{T}_{2}) \, \text{vec}(F^{T}) = 0.$$ Denote by $K$ the commutation matrix of $F$'s dimension. Then $$(P_{2}^{T} \otimes P_{1}^{T})\, \text{vec}(F) + (P_{1}^{T} \otimes P^{T}_{2}) \, K^{T} \, \text{vec}(F) = 0.$$

Finally, $F$ can be solved for by unvectorizing the solution to $$\left((P_{2}^{T} \otimes P_{1}^{T}) + (P_{1}^{T} \otimes P^{T}_{2}) \, K^{T} \right) \, \text{vec}(F) = 0,$$ which you can solve with scipy. Of course, the commutation matrix should be computed once and used throughout the for loop. I prefer this approach because it works with the camera matrices directly. I suspect it will be faster but cannot confirm.

Note also that the fundamental matrix between cameras $i$ and $j$ is the transpose of the fundamental matrix between cameras $j$ and $i$. This allows you to cut the number of fundamental matrices you are computing in half.