The `pfaffian` function¶

The Pfaffian of a \(2n \times 2n\) skew-symmetric matrix \(A\) (a square matrix with \(A^{\top} = -A\)) is a polynomial in its entries, with several equivalent definitions.

Formal definition (sum over perfect matchings). Let \(\Pi\) be the set of all perfect matchings of \(\{1, 2, \dots, 2n\}\): partitions into \(n\) unordered pairs \(\alpha = \{(i_1, j_1), \dots, (i_n, j_n)\}\) with every \(i_k < j_k\). Then

\[\operatorname{Pf}(A) = \sum_{\alpha \in \Pi} \operatorname{sgn}(\alpha)\, a_{i_1 j_1}\, a_{i_2 j_2} \cdots a_{i_n j_n},\]

where \(\operatorname{sgn}(\alpha)\) is the sign of the permutation that sends \((1, 2, \dots, 2n)\) to \((i_1, j_1, \dots, i_n, j_n)\). Equivalently, as a sum over the full symmetric group \(S_{2n}\),

\[\operatorname{Pf}(A) = \frac{1}{2^n\, n!} \sum_{\sigma \in S_{2n}} \operatorname{sgn}(\sigma)\, \prod_{k=1}^{n} a_{\sigma(2k-1)\, \sigma(2k)}.\]

Recursive definition (expansion along a row). The Pfaffian also satisfies a Laplace-like recursion. Writing \(A_{\hat{1}\hat{\jmath}}\) for the \((2n-2) \times (2n-2)\) matrix obtained by deleting rows and columns \(1\) and \(j\) from \(A\),

\[\operatorname{Pf}(A) = \sum_{j=2}^{2n} (-1)^{j}\, a_{1 j}\, \operatorname{Pf}\!\big(A_{\hat{1}\hat{\jmath}}\big),\]

with the base cases \(\operatorname{Pf} = 1\) for the empty \(0 \times 0\) matrix and \(\operatorname{Pf}\!\begin{pmatrix} 0 & a \\ -a & 0 \end{pmatrix} = a\) for the \(2 \times 2\) case.

Relation to the determinant. All of these agree, and the Pfaffian is the signed square root of the determinant:

\[\operatorname{Pf}(A)^2 = \det(A).\]

Where \(\sqrt{\det A}\) loses the sign, the Pfaffian keeps it. torch_pfaffian.pfaffian computes it for PyTorch tensors. This tutorial shows the inputs the function accepts and the outputs it returns.

import torch

from torch_pfaffian import pfaffian

# Seed the RNG so the random example matrices below are reproducible.
_ = torch.manual_seed(0)

Input: a single skew-symmetric matrix¶

The basic input is a skew-symmetric matrix of shape (2n, 2n), and the output is a scalar tensor: its Pfaffian. For a \(2 \times 2\) matrix \(\begin{pmatrix} 0 & a \\ -a & 0 \end{pmatrix}\) the Pfaffian is simply \(a\).

matrix = torch.tensor([[0.0, -3.0], [3.0, 0.0]])
pf = pfaffian(matrix)

print("input shape :", tuple(matrix.shape))
print("output      :", pf)

input shape : (2, 2)
output      : tensor(-3.)

Any skew-symmetric matrix works. A quick way to build one is A = M - M.T, which is skew-symmetric by construction. The defining identity \(\operatorname{Pf}(A)^2 = \det(A)\) then holds:

m = torch.randn(6, 6)
a = m - m.T  # skew-symmetric: a.T == -a

pf = pfaffian(a)
print("Pf(A)      =", pf.item())
print("Pf(A) ** 2 =", (pf**2).item())
print("det(A)     =", torch.linalg.det(a).item())

Pf(A)      = 0.682274341583252
Pf(A) ** 2 = 0.46549826860427856
det(A)     = 0.4654996395111084

Input: a batch of matrices¶

Leading dimensions are treated as a batch. An input of shape (..., 2n, 2n) returns one Pfaffian per matrix, with shape (...,).

batch = torch.randn(4, 8, 8)
batch = batch - batch.transpose(-1, -2)  # skew-symmetric in the last two dims

pfs = pfaffian(batch)
print("input shape :", tuple(batch.shape))
print("output shape:", tuple(pfs.shape))
print("output      :", pfs)

input shape : (4, 8, 8)
output shape: (4,)
output      : tensor([-22.6322,  11.7543,   5.5873, -16.7749])

Input dtype: real or complex¶

The Pfaffian is defined for complex skew-symmetric matrices too. The output matches the input’s dtype: a real input gives a real Pfaffian, a complex input gives a complex one.

real = torch.randn(4, 4)
real = real - real.T

complex_matrix = torch.randn(4, 4, dtype=torch.complex64)
complex_matrix = complex_matrix - complex_matrix.transpose(-1, -2)

print("real input    ->", pfaffian(real).dtype, "|", pfaffian(real).item())
print("complex input ->", pfaffian(complex_matrix).dtype, "|", pfaffian(complex_matrix).item())

real input    -> torch.float32 | -3.7162179946899414
complex input -> torch.complex64 | (-0.5062580108642578+0.7343237400054932j)

Output: signed value or magnitude¶

By default pfaffian returns the signed Pfaffian. Pass sign=False to get only its magnitude \(|\operatorname{Pf}(A)|\), which takes a cheaper determinant-based path when the sign is not needed. Here the \(2 \times 2\) matrix has a negative Pfaffian, so the two differ by the sign:

print("signed (default) :", pfaffian(matrix).item())
print("magnitude        :", pfaffian(matrix, sign=False).item())

signed (default) : -3.0
magnitude        : 3.0

Output: differentiable¶

If the input requires gradients, so does the output: pfaffian is fully compatible with torch.autograd, so it composes with any PyTorch computation and backpropagates.

a_grad = a.clone().requires_grad_(True)
pfaffian(a_grad).backward()

print("input.requires_grad :", a_grad.requires_grad)
print("gradient shape      :", tuple(a_grad.grad.shape))

input.requires_grad : True
gradient shape      : (6, 6)

Summary¶

Input to `pfaffian`	Output
skew-symmetric matrix `(2n, 2n)`	scalar Pfaffian
batch `(..., 2n, 2n)`	Pfaffians of shape `(...,)`
real or complex dtype	Pfaffian of the same dtype
`requires_grad=True`	differentiable output
CPU or CUDA tensor	result on the same device

In every case the output shares the input’s dtype, device, and backend, and returns a value whose square is \(\det(A)\).

The pfaffian function¶