torch_pfaffian.strategies package¶

Submodules¶

torch_pfaffian.strategies.pfaffian_block_det module¶

class torch_pfaffian.strategies.pfaffian_block_det.PfaffianBlockDet(*args, **kwargs)¶

Bases: PfaffianStrategy

Pfaffian of a block-antidiagonal skew-symmetric matrix from the determinant of its off-diagonal block.

For a matrix of the form

\[\begin{split}A = \begin{pmatrix} 0 & B \\ -B^T & 0 \end{pmatrix}, \qquad \text{pf}(A) = (-1)^{n(n-1)/2} \det(B),\end{split}\]

with B the upper-right n x n block, a single n x n determinant gives the signed Pfaffian. This identity holds only for that block form: on a general skew-symmetric matrix the formula returns a wrong value. The forward therefore validates the input with is_block_form() (both diagonal blocks exactly zero, even dimension) and raises ValueError otherwise. The check reads both diagonal blocks and costs up to three host synchronizations; set CHECK_BLOCK_FORM to False, or pass False as the second positional argument of apply() (apply(matrix, False)), to skip it for inputs known to have the block form.

The backward returns the same gradient convention as the general strategies, (1 / 2) pf(A) (A^{-1})^T on the full matrix, split between the upper-right and lower-left blocks.

CHECK_BLOCK_FORM = True¶

NAME = 'PfaffianBlockDet'¶

SVD_COFACTOR_MIN_BLOCK_SIZE = 17¶

static backward(ctx: BackwardCFunction, grad_output: Tensor) → tuple[torch.Tensor | None, None]¶

Compute the gradient of the Pfaffian with respect to the full input matrix.

With pf = c * det(B), c = (-1)^{n(n-1)/2} and B the upper-right block, Jacobi’s formula gives d pf / d B = c * adj(B)^T where adj(B) = det(B) B^{-1} is the classical adjugate. For the block form, A^{-1} = [[0, -B^{-T}], [B^{-1}, 0]], so the general-strategy gradient (1 / 2) pf(A) (A^{-1})^T is (1 / 2) c adj(B)^T in the upper-right block, its negated transpose in the lower-left block, and zero in the diagonal blocks. Returning that split (rather than the whole derivative in the upper-right block) makes A.grad independent of which strategy computed the Pfaffian; the gradient with respect to B of A = [[0, B], [-B^T, 0]] is the same either way.

For an invertible block the adjugate is the single inverse pf * (B^{-1})^T. For a (numerically) singular block that product is inaccurate or garbage, yet the true derivative c * adj(B)^T is finite and generally nonzero. The adjugate is then computed exactly by _exact_cofactor_matrix() (from an SVD, with the cofactor determinants as fallback), on the flagged batch elements only so invertible inputs keep the single cheap inverse. A block is flagged by PfaffianStrategy._flagged_inverse(), the same dtype-adaptive criteria as PfaffianStrategy.pfaffian_grad_matrix() (LU pivot ratio, LU info and inverse residual; an exactly singular block can have a round-off nonzero forward determinant, so pf == 0 alone misses it). Flagging a well-conditioned block only costs speed, never accuracy. Flagged blocks whose grad_output is exactly zero get a zero gradient without the exact path, since the gradient is linear in grad_output. With EXACT_SINGULAR_GRAD set to False the backward is host-synchronization-free and blocks singular at round-off level (see PfaffianStrategy._round_off_inverse()) get a zero gradient. That test is an LU pivot ratio below n * eps, so an invertible but strongly graded block (for example diag(1, 1e-9) in float32) is also given a zero gradient: this mode trades accuracy for speed.

The Pfaffian is holomorphic in the entries of A, so for complex inputs the backward returns the conjugate of the analytic derivative times grad_output, which is PyTorch’s Wirtinger convention. For real inputs the conjugation is a no-op.

Parameters:

• ctx – Context holding the saved input matrix and the forward Pfaffian.

• grad_output – Gradient of the output with respect to the loss.

Returns:

Gradient of the input matrix (or None when the input does not require grad), and None for the check_block_form flag.

Return type:

tuple[torch.Tensor | None, None]

static forward(matrix: Tensor, check_block_form: bool | None = None) → Tensor¶

Compute the Pfaffian of a block-antidiagonal skew-symmetric matrix from its upper-right block.

Parameters:

• matrix – Skew-symmetric matrix [[0, B], [-B^T, 0]] of shape (..., 2n, 2n).

• check_block_form – Whether to validate the block form. None (default) uses CHECK_BLOCK_FORM.

Returns:

The signed Pfaffian, of shape (...,).

Return type:

torch.Tensor

Raises:

ValueError – If the check is enabled and matrix does not have the block form.

static is_block_form(matrix: Tensor) → bool¶

Whether every matrix of the batch has the form [[0, B], [-B^T, 0]] with square B.

Only the structure is checked (square, even dimension, both diagonal blocks exactly zero); the skew-symmetry relation between the two off-diagonal blocks is assumed, as for every strategy. The entry (0, 1) is tested first so a general skew-symmetric input is rejected after reading a single entry per batch element.

Parameters:

matrix – Matrix of shape (..., m, m).

Returns:

True when all batch elements have the block-antidiagonal form.

Return type:

bool

static setup_context(ctx: FunctionCtx, inputs, output)¶

torch_pfaffian.strategies.pfaffian_det module¶

class torch_pfaffian.strategies.pfaffian_det.PfaffianDet(*args, **kwargs)¶

Bases: PfaffianStrategy

Compute the Pfaffian as pf = sqrt(|det(A)|) and rely on autograd for the backward pass.

The other strategies are full torch.autograd.Function implementations with a custom analytic backward. This strategy instead differentiates straight through det and sqrt with PyTorch autograd, so it overrides apply() to bypass the torch.autograd.Function machinery. The radicand is clamped to EPSILON for numerical stability, matching the other determinant-based strategy.

NAME = 'PfaffianDet'¶

static apply(matrix: Tensor) → Tensor¶

torch_pfaffian.strategies.pfaffian_fdbpf module¶

class torch_pfaffian.strategies.pfaffian_fdbpf.PfaffianFDBPf(*args, **kwargs)¶

Bases: PfaffianStrategy

This class implements the Pfaffian using the determinant of the matrix for the forward pass and the derivative of the Pfaffian with respect to the input matrix for the backward pass.

NAME = 'PfaffianFDBPf'¶

static backward(ctx: BackwardCFunction, grad_output: Tensor) → Tensor | None¶

Gradient of the Pfaffian magnitude with respect to the input matrix.

\[\frac{\partial |\text{pf}(A)|}{\partial A_{ij}} = \frac{|\text{pf}(A)|}{2} (A^{-1})_{ji}\]

For complex inputs the magnitude is a real function of a complex matrix, so the backward follows PyTorch’s Wirtinger convention (z.grad = 2 d L / d conj(z) for a real loss) and returns (|pf(A)| / 2) conj(A^{-T}). For real inputs the conjugation is a no-op.

The inverse uses torch.linalg.inv() (an LU factorization) rather than torch.linalg.pinv() (an SVD), which can fail to converge on ill-conditioned or near-repeated-singular-value inputs. The forward floors the magnitude at sqrt(EPSILON), so singular elements (pf at that floor, and the always-singular odd-dimensional inputs whose pf is 0) are replaced by the identity before the batched inverse: inv raises on an exactly-singular matrix, and there the gradient is negligible anyway since it scales with pf.

Parameters:

• ctx – Context holding the saved input matrix and the forward Pfaffian magnitude.

• grad_output – Gradient of the output with respect to the loss.

Returns:

Gradient of the input matrix, or None when the input does not require grad.

Return type:

torch.Tensor | None

static forward(matrix: Tensor)¶

torch_pfaffian.strategies.pfaffian_masked_magnitude module¶

class torch_pfaffian.strategies.pfaffian_masked_magnitude.PfaffianMaskedMagnitude(*args, **kwargs)¶

Bases: PfaffianStrategy

Pfaffian magnitude sqrt(max(|det(A)|, epsilon)) whose gradient is zero where Pf(A) vanishes up to round-off.

The forward is the fast determinant path of PfaffianDet, floored at sqrt(epsilon) with a per-call epsilon. The backward returns (|Pf(A)| / 2) conj(A^{-T}) (PyTorch’s Wirtinger convention for a real output of a complex input), except that it gives a zero gradient to:

• elements whose determinant was clamped to epsilon in the forward (the mask is saved in the forward rather than recomputed from the floored value, which could be off by one ulp);

• elements whose LU factorization fails, or whose 1-norm reciprocal condition number is not finite;

• elements that vanish up to round-off, flagged when both rcond_1(A) <= ROUNDOFF_RCOND_SCALE * eps and sigma_min(A) / sigma_max(A) <= ROUNDOFF_SPECTRAL_SCALE * eps, with eps the machine epsilon of the dtype.

A Pfaffian that is zero in exact arithmetic is computed as round-off, and the unmasked backward then divides noise by noise and returns an O(1) gradient. No numerical test separates such a structural zero from a genuine value at round-off level (in float32 the rare outcomes of a state near a basis state are genuine yet already 10 to 35% off), so this strategy is opt-in: use it when exact zeros are expected, for example the parity-forbidden outcomes of a basis state.

The determinant is computed with a backward error ~ eps * sigma_max, so the round-off of |Pf(A)| = prod sigma_pair is ~ eps * sigma_max * |Pf(A)| / sigma_min: the Pfaffian is unresolved when sigma_min / sigma_max is at round-off level, at every size. For a skew matrix rcond_1(A) <= sigma_min / sigma_max, and the 1-norm ratio can be up to m times smaller, so the 1-norm test alone would mask small but resolved Pfaffians, while the spectral test alone would miss the structural zeros of long circuits. Both tests must agree:

• a structural zero rounded once measures sigma_min / sigma_max <= 0.5 eps up to m = 256; accumulated round-off (up to 1e5 Givens rotations applied in place in complex128) raises it to ~10 eps while its rcond_1 stays below 4 eps. Longer circuits drift past both tests;

• a Pfaffian with sigma_min / sigma_max ~ 4 eps is still accurate to a few percent, so only the 1-norm test separates it from a structural zero.

The backward costs one inverse, two 1-norms and one set of singular values per element. When the SVD fails to converge, the whole batch falls back on the 1-norm test alone.

The criterion has a known limit: it is not scale-free across blocks, so a resolved element whose exact block structure spans a very wide range (about 0.25 / eps between block scales) can be masked; once rotated, such an element is no longer resolved in floating point. Elements whose determinant overflows are recovered by torch_pfaffian.pfaffian() in the log domain, which carries the unmasked gradient. Only dtypes with an LU factorization are supported (float32, float64, complex64, complex128).

Odd-dimensional matrices are singular, so their magnitude is the floor sqrt(epsilon) with a zero gradient. apply() returns only the magnitude.

NAME = 'PfaffianMaskedMagnitude'¶

ROUNDOFF_RCOND_SCALE = 4.0¶

ROUNDOFF_SPECTRAL_SCALE = 16.0¶

classmethod apply(matrix: Tensor, epsilon: float | None = None) → Tensor¶

Compute the masked-gradient Pfaffian magnitude.

Parameters:

• matrix – Skew-symmetric matrices of shape (..., m, m), real or complex.

• epsilon – Floor of |det(A)|. None (default) uses PfaffianStrategy.EPSILON.

Returns:

The magnitude of shape (...,) in the real dtype of matrix.

Return type:

torch.Tensor

static backward(ctx: BackwardCFunction, grad_output: Tensor, grad_clamped: Tensor | None = None) → tuple[torch.Tensor | None, None]¶

Gradient grad_output * (|Pf(A)| / 2) conj(A^{-T}), zero at clamped, singular and round-off elements.

Parameters:

• ctx – Context holding the input matrix, the forward magnitude and the clamped mask.

• grad_output – Gradient of the loss with respect to the magnitude, of shape (...,).

• grad_clamped – Unused; the clamped mask is not differentiable.

Returns:

Gradient of the input matrix (or None when it does not require grad), and None for epsilon.

Return type:

tuple[torch.Tensor | None, None]

static forward(matrix: Tensor, epsilon: float | None = None) → tuple[torch.Tensor, torch.Tensor]¶

Compute the floored magnitude and the mask of the elements clamped to the floor.

Parameters:

• matrix – Skew-symmetric matrices of shape (..., m, m), real or complex.

• epsilon – Floor of |det(A)|. None (default) uses PfaffianStrategy.EPSILON.

Returns:

The magnitude sqrt(max(|det(A)|, epsilon)) of shape (...,) in the real dtype of matrix, and the boolean clamped mask of shape (...,).

Return type:

tuple[torch.Tensor, torch.Tensor]

static setup_context(ctx: FunctionCtx, inputs, output)¶

torch_pfaffian.strategies.pfaffian_parlett_reid module¶

class torch_pfaffian.strategies.pfaffian_parlett_reid.PfaffianParlettReid(*args, **kwargs)¶

Bases: PfaffianStrategy

Compute the signed Pfaffian of a skew-symmetric matrix with the Parlett-Reid algorithm.

The forward performs a batched skew-tridiagonalization with partial pivoting: the batch is processed together and only the n / 2 column-elimination steps are sequential. The result is the signed Pfaffian, unlike the determinant-based strategies which return only the magnitude. The backward uses the closed form d pf(A) / d A = (1 / 2) pf(A) (A^{-1})^T via the Pfaffian adjugate, with no autograd graph over the elimination. For invertible inputs this is a single inverse; for singular inputs (pf == 0) the adjugate is computed exactly from minor Pfaffians, so the gradient is correct everywhere (see PfaffianStrategy.pfaffian_grad_matrix()).

The input is a skew-symmetric matrix of shape (..., 2n, 2n).

NAME = 'PfaffianParlettReid'¶

static backward(ctx: BackwardCFunction, grad_output: Tensor) → Tensor | None¶

Gradient of the signed Pfaffian with respect to the input matrix.

\[\frac{\partial \text{pf}(A)}{\partial A_{ij}} = \frac{\text{pf}(A)}{2} (A^{-1})_{ji}\]

Parameters:

• ctx – Context holding the saved input matrix and the forward Pfaffian.

• grad_output – Gradient of the output with respect to the loss.

Returns:

Gradient of the input matrix, or None when the input does not require grad.

Return type:

torch.Tensor | None

static forward(matrix: Tensor) → Tensor¶

torch_pfaffian.strategies.pfaffian_rust_parlett_reid module¶

class torch_pfaffian.strategies.pfaffian_rust_parlett_reid.RustPfaffianParlettReid(*args, **kwargs)¶

Bases: PfaffianStrategy

Compute the signed Pfaffian with the Rust Parlett-Reid kernel.

The forward moves the input to a contiguous CPU array and calls the compiled torch_pfaffian._rust kernel for the input dtype, each running natively at that precision: float16, float32/complex64 and float64/complex128 map to the half-, single- and double-precision kernels. Half precision is the caller’s explicit choice and carries its risk: the elimination runs entirely in float16, which is the least accurate kernel and can overflow the narrow float16 range to inf for larger or larger-scaled matrices (use float32 for a measurably more accurate result). Complex inputs are computed natively over the complex field (no imaginary part is discarded), giving the correct complex signed Pfaffian. Any other dtype raises TypeError. The result is cast back to the input dtype and device. The backward is the same as PfaffianParlettReid (the Pfaffian adjugate d pf(A) / d A = (1 / 2) pf(A) (A^{-1})^T, exact for invertible and singular inputs and conjugated for complex autograd), computed in PyTorch. CUDA inputs are evaluated on CPU for the forward; the backward runs on the input device.

The input is a skew-symmetric matrix of shape (..., n, n).

NAME = 'RustPfaffianParlettReid'¶

static backward(ctx: BackwardCFunction, grad_output: Tensor) → Tensor | None¶

Gradient of the signed Pfaffian with respect to the input matrix.

\[\frac{\partial \text{pf}(A)}{\partial A_{ij}} = \frac{\text{pf}(A)}{2} (A^{-1})_{ji}\]

Parameters:

• ctx – Context holding the saved input matrix and the forward Pfaffian.

• grad_output – Gradient of the output with respect to the loss.

Returns:

Gradient of the input matrix, or None when the input does not require grad.

Return type:

torch.Tensor | None

static forward(matrix: Tensor) → Tensor¶

torch_pfaffian.strategies.pfaffian_slog module¶

class torch_pfaffian.strategies.pfaffian_slog.SlogPfaffianStrategy(*args, **kwargs)¶

Bases: Function

Log-domain signed Pfaffian returning (phase, log|Pf|) with an analytic backward for log|Pf|.

This is a standalone torch.autograd.Function (deliberately not a PfaffianStrategy subclass, so it does not enter the scalar-returning strategy registry). The forward is the host-synchronization-free parlett_reid_slog(); the returned phase is marked non-differentiable (for real inputs it is piecewise constant, and callers needing the full complex derivative use the linear-domain PfaffianParlettReid). log|Pf| is differentiable wherever Pf != 0, with gradient 0.5 * (A^{-1})^T under PyTorch’s Wirtinger convention (the .conj() is a no-op for real inputs); batch elements with Pf = 0 receive an exactly zero gradient.

The input is a skew-symmetric matrix of shape (..., m, m).

static backward(ctx: Any, *grad_outputs: Any) → Tensor¶

Gradient of log|Pf(A)| with respect to the input matrix.

\[\frac{\partial \log|\text{pf}(A)|}{\partial A_{ij}} = \frac{1}{2} (A^{-1})_{ji}\]

The gradient of the non-differentiable phase output is ignored. Singular elements (phase == 0) are replaced by the identity before the batched inverse and then assigned an exactly zero gradient, so the backward never synchronizes with the host and never raises.

Parameters:

• ctx – Context holding the saved input matrix and the forward phase.

• grad_outputs – Gradients of (phase, log_abs) with respect to the loss; only the log_abs gradient is used.

Returns:

Gradient of the input matrix, of shape (..., m, m).

Return type:

torch.Tensor

static forward(matrix: Tensor) → tuple[torch.Tensor, torch.Tensor]¶

static setup_context(ctx: Any, inputs: tuple[torch.Tensor, ...], output: tuple[torch.Tensor, torch.Tensor]) → None¶

torch_pfaffian.strategies.pfaffian_small module¶

class torch_pfaffian.strategies.pfaffian_small.PfaffianSmall(*args, **kwargs)¶

Bases: PfaffianStrategy

Signed Pfaffian for small even dimensions via the unrolled perfect-matching kernel.

Like PfaffianDet, this strategy differentiates straight through plain tensor ops rather than a custom analytic backward, so it overrides apply() to call small_pfaffian() directly and lets autograd handle the backward. The kernel is exact and autograd-friendly for even dimensions m <= SMALL_N_MAX.

For m >= 2 the input is re-skew-symmetrized (A = 0.5 * (A - A^T)) before the kernel. The unrolled kernel only reads the upper triangle, so raw autograd would pile the whole gradient there, whereas the analytic vector-Jacobian product of the other strategies returns the skew-split convention 0.5 * pf(A) (A^{-1})^T. The projection is idempotent for skew inputs, so the forward is unchanged while the gradient matches the other strategies (including at m = 2, where the gradient is [[0, 0.5], [-0.5, 0]] rather than the upper-triangle [[0, 1], [0, 0]]).

The input is a skew-symmetric matrix of shape (..., m, m).

NAME = 'PfaffianSmall'¶

static apply(matrix: Tensor) → Tensor¶

torch_pfaffian.strategies.strategy module¶

class torch_pfaffian.strategies.strategy.PfaffianStrategy(*args, **kwargs)¶

Bases: Function

EPSILON = 1e-30¶

EXACT_SINGULAR_GRAD = True¶

INVERSE_RESIDUAL_RTOL_EXPONENT = 0.5¶

MINOR_CHUNK_ELEMENTS = 1048576¶

NAME = 'PfaffianStrategy'¶

NULL_SPACE_SEPARATION_EXPONENT = 0.25¶

SINGULARITY_RTOL_EXPONENT = 0.75¶

static backward(ctx: BackwardCFunction, grad_output)¶

static forward(matrix: Tensor)¶

classmethod pfaffian_grad_matrix(matrix: Tensor, pfaffian: Tensor, grad_output: Tensor) → Tensor¶

Gradient of the signed Pfaffian with respect to the input matrix.

Uses the closed form d pf(A) / d A = (1 / 2) pf(A) (A^{-1})^T via the Pfaffian adjugate pf(A) A^{-1}. For well-conditioned inputs the adjugate is pf(A) * inv(A) (a single inverse); for (numerically) singular inputs, where that product would be inaccurate or garbage, the adjugate is recomputed exactly by _singular_pfaffian_adjugate() (from the singular values and at most one minor Pfaffian when the corank is clear, else from all the minor Pfaffians of _pfaffian_adjugate(), using cls’s own forward). The minor-based path runs only on the flagged batch elements, so well-conditioned inputs keep the single cheap inverse. Flagged elements whose grad_output is exactly zero (for example the elements of a magnitude whose gradient vanishes at pf = 0) get a zero gradient without the minor-based path, since the gradient is linear in grad_output; this skip is disabled when grad_output itself requires grad (double backward, forward-mode products such as torch.autograd.functional.jvp).

An element is flagged by _flagged_inverse() (the LU pivot ratio min|U_ii| / max|U_ii| <= eps^0.75, a nonzero LU info, or an inverse residual ||A A^{-1} - I||_max > eps^0.5), or when its forward Pfaffian is exactly 0 (which also catches a Pfaffian that underflowed on a well-conditioned matrix with tiny entries). An exactly-singular matrix has a round-off-sized LU pivot (and a round-off forward Pfaffian, rarely exactly 0), so the pivot ratio catches it, while the residual catches ill-conditioned elements whose pivots look benign. The pivot ratio is scale-free: unlike a test of |pf| against max|A_ij|^(n/2), it does not flag well-conditioned matrices whose entries span several orders of magnitude (e.g. block-diagonal matrices). Flagging a well-conditioned element only costs speed, never accuracy.

The inverse is an LU factorization (torch.linalg.lu_factor_ex() and torch.linalg.lu_solve(), or torch.linalg.inv_ex() in double backward) rather than torch.linalg.pinv() (an SVD), which can fail to converge on ill-conditioned or near-repeated-singular-value inputs. A skew-symmetric A is invertible exactly when pf(A) != 0 (since det(A) = pf(A)^2), so the inverse is only relied upon on the unflagged elements. Flagged elements are solved against the identity factors instead, so the batched solve never divides by a zero pivot.

The Pfaffian is holomorphic in the entries of A, so for complex inputs the backward returns the conjugate of the analytic derivative, conj(d pf / d A) * grad_output, which is PyTorch’s Wirtinger convention for complex autograd. For real inputs the conjugation is a no-op, so real gradients are unchanged.

All of the above assumes EXACT_SINGULAR_GRAD is True (the default). When it is False the backward takes a fully host-synchronization-free path, skipping the .any() branch and the adjugate entirely: elements singular at round-off level (an LU pivot ratio at most n * eps, an LU failure, or pf == 0; finite inputs only, so NaN still propagates) receive an exactly zero gradient (the true derivative there is nonzero only at corank exactly 2), and every other element, however ill-conditioned, the single inverse pf(A) A^{-1} (accurate to about cond(A) * eps). In this mode flagging trades accuracy for speed, unlike the exact mode.

Parameters:

• matrix – The saved input matrix of shape (..., n, n).

• pfaffian – The saved forward Pfaffian of shape (...,).

• grad_output – Gradient of the output with respect to the loss, of shape (...,).

Returns:

Gradient of the input matrix, of shape (..., n, n).

Return type:

torch.Tensor

static setup_context(ctx: FunctionCtx, inputs, output)¶

Module contents¶