torch_pfaffian.low_rank package¶

Submodules¶

torch_pfaffian.low_rank.low_rank_pfaffian module¶

class torch_pfaffian.low_rank.low_rank_pfaffian.LowRankPfaffian¶

Bases: object

Ratios of Pfaffians and inverse updates under skew-symmetric low-rank updates.

For a skew-symmetric invertible A and u of shape (n, 2k), the updated matrix is

\[\begin{split}A' = A - u J u^T, \qquad J = \begin{pmatrix} 0 & I_k \\ -I_k & 0 \end{pmatrix}\end{split}\]

(the convention of lrux, arXiv:2602.05255), and with the 2k x 2k skew-symmetric matrix R = J + u^T A^{-1} u,

\[\frac{\mathrm{Pf}(A')}{\mathrm{Pf}(A)} = \frac{\mathrm{Pf}(R)}{\mathrm{Pf}(J)}, \qquad A'^{-1} = A^{-1} + (A^{-1} u) R^{-1} (A^{-1} u)^T.\]

The ratio costs O(n^2 k) instead of the O(n^3) of a new Pfaffian, and so does the updated inverse (Woodbury identity), which is skew-symmetric up to round-off. A rank-2 update u = (x, y) adds y x^T - x y^T to A; replacing one row and column (ratio_replace()) costs O(n) for the ratio.

Only tensor operations are used, so every method is batched, runs on any device without host synchronization, and is differentiable. Complex inputs use the transpose, not the conjugate transpose, as the Pfaffian is bilinear. float16, bfloat16 and complex32 inputs are computed in float32 or complex64 with a UserWarning, as in pfaffian(). When a ratio is zero, the updated matrix and R are singular: the returned inverse is then meaningless (huge or non-finite) and should be recomputed from the updated matrix when it is needed.

Gradients are exact for any skew-symmetric construction of the inputs (e.g. from circuit parameters). The gradient with respect to a raw inverse leaf is projected onto skew-symmetric matrices (SkewGradient), the convention of pfaffian(). The gradients with respect to update, change, border and corner sit on the entries that are read: when they are slices of a larger skew-symmetric leaf matrix, only the skew-symmetric part of that matrix’s gradient is meaningful.

PfaffianUpdater chains such updates, with delayed inverse updates.

static add_product(matrix: Tensor, left: Tensor, right: Tensor) → Tensor¶

matrix + left @ right in one fused pass over matrix when the batch shapes allow it.

The low-rank corrections of the inverse are skew-symmetric by construction, so the result stays skew-symmetric up to round-off without an explicit (and, at large n, costly) symmetrization.

Parameters:

• matrix – Matrix of shape (..., n, n).

• left – Left factor, of shape (..., n, p).

• right – Right factor, of shape (..., p, n).

Returns:

matrix + left @ right, of shape (..., n, n).

Return type:

torch.Tensor

static as_index(index: Tensor | int, device: device) → Tensor¶

Row index of ratio_replace() as an int64 tensor on device.

Parameters:

• index – An integer or an integer tensor.

• device – Device of the matrices.

Returns:

The index tensor.

Return type:

torch.Tensor

Raises:

TypeError – If index is not an integer.

static project(inverse: Tensor, vectors: Tensor) → Tensor¶

Product A^{-1} v with a skew-symmetric A^{-1}, computed as -(v^T A^{-1})^T.

The transposed product reads A^{-1} row by row and is about 3x faster on CPU for the skinny v of low-rank updates; the two agree up to the round-off in the skew-symmetry of A^{-1}.

Parameters:

• inverse – Skew-symmetric matrix A^{-1}, of shape (..., n, n).

• vectors – Vectors v, of shape (..., n, p).

Returns:

A^{-1} v, of shape (..., n, p).

Return type:

torch.Tensor

static ratio(inverse: Tensor, update: Tensor) → Tensor¶

Ratio Pf(A - u J u^T) / Pf(A) from the inverse A^{-1}.

Parameters:

• inverse – Inverse A^{-1} of the original skew-symmetric matrix, of shape (..., n, n).

• update – Update vectors u, of shape (..., n, 2k); the batch dimensions broadcast with those of inverse.

Returns:

The ratio, of shape (...).

Return type:

torch.Tensor

Raises:

• ValueError – If the shapes do not match or update has an odd or zero number of columns.

• TypeError – If the dtypes or devices of the inputs differ, or the dtype is not supported.

static ratio_and_inverse(inverse: Tensor, update: Tensor) → tuple[torch.Tensor, torch.Tensor]¶

Ratio Pf(A - u J u^T) / Pf(A) of ratio() and the inverse (A - u J u^T)^{-1}.

Parameters:

• inverse – Inverse A^{-1} of the original skew-symmetric matrix, of shape (..., n, n).

• update – Update vectors u, of shape (..., n, 2k); the batch dimensions broadcast with those of inverse.

Returns:

The ratio, of shape (...), and the updated inverse, of shape (..., n, n).

Return type:

tuple[torch.Tensor, torch.Tensor]

Raises:

• ValueError – If the shapes do not match or update has an odd or zero number of columns.

• TypeError – If the dtypes or devices of the inputs differ, or the dtype is not supported.

static ratio_and_r_inverse(update: Tensor, projected: Tensor) → tuple[torch.Tensor, torch.Tensor]¶

Ratio of ratio() and the inverse of its matrix R, from the projection A^{-1} u.

Parameters:

• update – Update vectors u, of shape (..., n, 2k).

• projected – Projection A^{-1} u, of shape (..., n, 2k).

Returns:

The ratio, of shape (...), and R^{-1}, of shape (..., 2k, 2k) (huge or non-finite when the ratio is zero).

Return type:

tuple[torch.Tensor, torch.Tensor]

static ratio_extend(inverse: Tensor, border: Tensor, corner: Tensor) → Tensor¶

Ratio of Pfaffians when a skew-symmetric matrix grows by 2k rows and columns.

For the bordered matrix

\[\begin{split}M = \begin{pmatrix} A & B \\ -B^T & C \end{pmatrix}, \qquad \frac{\mathrm{Pf}(M)}{\mathrm{Pf}(A)} = \mathrm{Pf}\left(C + B^T A^{-1} B\right),\end{split}\]

the Pfaffian of the 2k x 2k Schur complement, in O(m^2 k) from A^{-1}. This is the step of sequential (qubit-by-qubit) sampling of matchgate circuits, where a qubit adds the 2 x 2 block of its two Majorana modes: the conditional probability of an outcome is |ratio| / 2. A may be empty (m = 0), when the ratio is Pf(C).

Only C depends on the outcome, so several candidate corners are evaluated at the cost of one: give them a leading dimension, e.g. corner of shape (2, ..., 2k, 2k) for the outcomes 0 and 1, and the ratios have shape (2, ...) (A^{-1} B is computed once, at the batch shape of inverse and border). PfaffianChain keeps A^{-1} from one step to the next.

Parameters:

• inverse – Inverse A^{-1} of the current skew-symmetric matrix, of shape (..., m, m).

• border – Coupling B of the current rows to the new ones, of shape (..., m, 2k).

• corner – Skew-symmetric block C of the new rows and columns, of shape (..., 2k, 2k).

Returns:

The ratio Pf(M) / Pf(A), of the broadcast batch shape (...).

Return type:

torch.Tensor

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices of the inputs differ, or the dtype is not supported.

static ratio_extend_and_inverse(inverse: Tensor, border: Tensor, corner: Tensor) → tuple[torch.Tensor, torch.Tensor]¶

Ratio of ratio_extend() and the inverse of the bordered matrix, in O(m^2 k).

With P = A^{-1} B and the Schur complement S = C + B^T A^{-1} B,

\[\begin{split}M^{-1} = \begin{pmatrix} A^{-1} + P S^{-1} P^T & -P S^{-1} \\ -S^{-1} P^T & S^{-1} \end{pmatrix}.\end{split}\]

Parameters:

• inverse – Inverse A^{-1} of the current skew-symmetric matrix, of shape (..., m, m).

• border – Coupling B of the current rows to the new ones, of shape (..., m, 2k).

• corner – Skew-symmetric block C of the new rows and columns, of shape (..., 2k, 2k).

Returns:

The ratio, of the broadcast batch shape (...), and M^{-1}, of shape (..., m + 2k, m + 2k) (huge or non-finite when the ratio is zero).

Return type:

tuple[torch.Tensor, torch.Tensor]

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices of the inputs differ, or the dtype is not supported.

static ratio_replace(inverse: Tensor, index: Tensor | int, change: Tensor) → Tensor¶

Ratio of Pfaffians when one row and the matching column of a skew-symmetric matrix change.

Row i of A becomes A[i] + c and column i becomes A[:, i] - c, so that A' stays skew-symmetric; c_i is ignored (the diagonal stays zero). This is the rank-2 update u = (c, e_i) of ratio(), with a one-hot e_i, for which

\[r = \frac{\mathrm{Pf}(A')}{\mathrm{Pf}(A)} = 1 + c^T A^{-1} e_i.\]

It reads one column of the inverse and costs O(n). To replace row i by a new row x (with x_i = 0), pass change = x - A[i]. ratio_replace_and_inverse() also updates the inverse.

Parameters:

• inverse – Inverse A^{-1} of the original skew-symmetric matrix, of shape (..., n, n).

• index – Index i of the row and column that change, an integer or an integer tensor of shape (...).

• change – Change c of row i, of shape (..., n).

Returns:

The ratio, whose shape is the broadcast of the batch shapes of inverse, index and change.

Return type:

torch.Tensor

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices of the inputs differ, index is not an integer tensor, or the dtype is not supported.

static ratio_replace_and_inverse(inverse: Tensor, index: Tensor | int, change: Tensor) → tuple[torch.Tensor, torch.Tensor]¶

Ratio of ratio_replace() and the inverse of the updated matrix, in O(n^2):

\[A'^{-1} = A^{-1} + \frac{(A^{-1} e_i)(A^{-1} c)^T - (A^{-1} c)(A^{-1} e_i)^T}{r}.\]

Parameters:

• inverse – Inverse A^{-1} of the original skew-symmetric matrix, of shape (..., n, n).

• index – Index i of the row and column that change, an integer or an integer tensor of shape (...).

• change – Change c of row i, of shape (..., n).

Returns:

The ratio, of the broadcast batch shape (...), and the updated inverse, of shape (..., n, n).

Return type:

tuple[torch.Tensor, torch.Tensor]

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices of the inputs differ, index is not an integer tensor, or the dtype is not supported.

static skew_identity(rank: int, *, dtype: dtype | None = None, device: device | str | None = None) → Tensor¶

Skew-symmetric identity J = [[0, I_k], [-I_k, 0]] of the low-rank Pfaffian updates.

Its Pfaffian is (-1)^(k (k - 1) / 2) (see skew_identity_pfaffian()) and its inverse is -J.

Parameters:

• rank – Half the size k of the matrix, at least 1.

• dtype – Dtype of the result (PyTorch’s default dtype when None).

• device – Device of the result (the default device when None).

Returns:

The matrix J, of shape (2k, 2k).

Return type:

torch.Tensor

Raises:

ValueError – If rank is smaller than 1.

static skew_identity_pfaffian(rank: int) → int¶

Pfaffian of skew_identity(), (-1)^(k (k - 1) / 2).

Parameters:

rank – Half the size k of the matrix.

Returns:

1 or -1.

Return type:

int

torch_pfaffian.low_rank.pfaffian_chain module¶

class torch_pfaffian.low_rank.pfaffian_chain.PfaffianChain(batch_shape: tuple[int, ...] | Size = (), *, dtype: dtype = torch.float64, device: device | str | None = None)¶

Bases: object

Pfaffian of a skew-symmetric matrix that grows block by block, with the ratio of each growth step.

The state is the inverse of the current m x m matrix A and its Pfaffian. A step appends 2k rows and columns, the coupling B (m x 2k) and the block C (2k x 2k):

\[\begin{split}\begin{pmatrix} A & B \\ -B^T & C \end{pmatrix}, \qquad \text{ratio} = \mathrm{Pf}\left(C + B^T A^{-1} B\right),\end{split}\]

in O(m^2 k) (see ratio_extend()), so the Pfaffians of all the leading blocks of an n x n matrix cost O(n^3) in total, against O(n^4) when each is computed anew.

This is the loop of qubit-by-qubit sampling of a matchgate circuit (as in MatchCake). The probability of the outcomes y_0 ... y_j is 2^{-(j+1)} |Pf(M_j)|, with M_j the principal submatrix of Lambda + Lambda_y on the Majorana modes of qubits 0 ... j. Qubit j appends its two modes: B is the covariance between the modes of the qubits sampled so far and its own, and C its 2 x 2 block, which depends on the outcome y_j. ratios() gives the ratio for every candidate block (the conditional probabilities are |ratio| / 2), and extend() appends the sampled one:

chain = PfaffianChain((shots,), dtype=covariance.dtype)
for qubit in range(n_qubits):
    modes = slice(2 * qubit, 2 * qubit + 2)
    border = covariance[..., : 2 * qubit, modes]                   # (shots, 2 qubit, 2)
    corners = covariance[..., modes, modes] + outcome_blocks       # (shots, 2, 2, 2), y = 0 and 1
    probabilities = chain.ratios(border, corners).abs() / 2        # (shots, 2)
    outcome = torch.multinomial(probabilities, 1)[..., 0]          # (shots,)
    chain.extend(border, corners[torch.arange(shots), outcome])

Every operation is a batched tensor operation on any device, without host synchronization, and is differentiable: gradients are exact for any skew-symmetric construction of the matrix (e.g. Lambda(t) from circuit parameters). The chain reads the blocks B and C and never -B^T, so when they are slices of a skew-symmetric leaf matrix, its gradient sits on those blocks and only its skew-symmetric part, g - g^T, matches that of pfaffian(). A zero ratio (an outcome of probability zero) makes the grown matrix singular: such a block must not be appended.

The batch dimensions of the borders and corners broadcast with batch_shape. The dtype is float32, float64, complex64 or complex128.

__init__(batch_shape: tuple[int, ...] | Size = (), *, dtype: dtype = torch.float64, device: device | str | None = None) → None¶

Start from the empty matrix (m = 0, Pfaffian 1).

Parameters:

• batch_shape – Batch shape of the chain, e.g. (shots,); the first step may broadcast it.

• dtype – Dtype of the matrices: float32, float64, complex64 or complex128.

• device – Device of the matrices.

Raises:

TypeError – If dtype is not supported.

property dimension: int¶

Size m of the current matrix.

Returns:

The number of rows of the current matrix.

Return type:

int

extend(border: Tensor, corner: Tensor) → Tensor¶

Append the rows and columns (B, C) to the current matrix.

Parameters:

• border – Coupling B of the current rows to the new ones, of shape (..., m, 2k).

• corner – Block C of the new rows and columns, of shape (..., 2k, 2k).

Returns:

The ratio of the step, Pf(M') / Pf(M), of shape (...).

Return type:

torch.Tensor

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices differ from those of the chain.

classmethod from_matrix(matrix: Tensor) → PfaffianChain¶

Start from a non-empty skew-symmetric matrix, computing its inverse and Pfaffian.

Parameters:

matrix – Invertible skew-symmetric matrix of shape (..., m, m), in float32, float64, complex64 or complex128.

Returns:

The chain, whose current matrix is matrix.

Return type:

PfaffianChain

Raises:

• ValueError – If matrix is not square.

• TypeError – If the dtype of matrix is not supported.

property inverse: Tensor¶

Inverse of the current matrix.

Returns:

The inverse, of shape (..., m, m).

Return type:

torch.Tensor

property pfaffian: Tensor¶

Pfaffian of the current matrix, the product of the ratios of the steps (times that of the initial matrix for from_matrix()).

Returns:

The Pfaffian, of shape (...).

Return type:

torch.Tensor

ratio(border: Tensor, corner: Tensor) → Tensor¶

Ratio Pf(M') / Pf(M) for one block, without changing the state.

Parameters:

• border – Coupling B of the current rows to the new ones, of shape (..., m, 2k).

• corner – Block C of the new rows and columns, of shape (..., 2k, 2k).

Returns:

The ratio, of shape (...).

Return type:

torch.Tensor

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices differ from those of the chain.

ratios(border: Tensor, corners: Tensor) → Tensor¶

Ratios Pf(M') / Pf(M) for several candidate blocks, without changing the state.

Parameters:

• border – Coupling B of the current rows to the new ones, of shape (..., m, 2k).

• corners – Candidate blocks C, of shape (..., c, 2k, 2k).

Returns:

The ratio for every candidate, of shape (..., c).

Return type:

torch.Tensor

Raises:

• ValueError – If the shapes do not match.

• TypeError – If the dtypes or devices differ from those of the chain.

torch_pfaffian.low_rank.pfaffian_updater module¶

class torch_pfaffian.low_rank.pfaffian_updater.PfaffianUpdater(matrix: Tensor, *, max_delay: int = 8)¶

Bases: object

Successive low-rank updates of a skew-symmetric matrix, with Pfaffian ratios and delayed inverse updates.

The state is the inverse A_0^{-1} of the last merged matrix and, for each of the tau updates accepted since, the projection a_t = A_{t-1}^{-1} u_t and R_t^{-1} (see LowRankPfaffian for the convention A_t = A_{t-1} - u_t J u_t^T). The current inverse is never formed between merges:

\[A_\tau^{-1} u = A_0^{-1} u + \sum_{t=1}^{\tau} a_t R_t^{-1} (a_t^T u),\]

so an update costs O(n^2 k + \tau n k^2) and does not rewrite the n x n inverse. After max_delay accepted updates the corrections are merged into the inverse. lrux (arXiv:2602.05255) suggests max_delay of about n / (20 k).

Every operation is a tensor operation on the batch, without host synchronization, so a batch of Monte Carlo walkers is updated at once: accept() takes a mask of the walkers whose move is accepted, and leaves the others unchanged. Gradients flow through the state like through any tensor computation.

Round-off accumulates over many updates, mostly in single precision. reset() recomputes the inverse from the current matrix (which the caller keeps); the running pfaffian is the Pfaffian of the initial matrix times the accepted ratios, and can over- or underflow over long runs.

The input matrix is a skew-symmetric matrix of shape (..., n, n), in float32, float64, complex64 or complex128.

__init__(matrix: Tensor, *, max_delay: int = 8) → None¶

Start from matrix, computing its inverse and its Pfaffian.

Parameters:

• matrix – Invertible skew-symmetric matrix of shape (..., n, n).

• max_delay – Number of accepted updates after which they are merged into the inverse, at least 1.

Raises:

• ValueError – If max_delay is smaller than 1 or matrix is not square.

• TypeError – If the dtype of matrix is not supported.

accept(update: Tensor, mask: Tensor | None = None) → Tensor¶

Apply the update A <- A - u J u^T to the batch elements selected by mask.

Parameters:

• update – Update vectors u, of shape (..., n, 2k).

• mask – Boolean tensor of shape (...), True where the update is applied; every element when None.

Returns:

The ratio of the applied update, of shape (...), equal to 1 where mask is False.

Return type:

torch.Tensor

accept_replace(index: Tensor | int, change: Tensor, mask: Tensor | None = None) → Tensor¶

Add change to row index and subtract it from column index, where mask is True.

Parameters:

• index – Index of the row and column, an integer or an integer tensor of shape (...).

• change – Change of the row, of shape (..., n).

• mask – Boolean tensor of shape (...), True where the change is applied; every element when None.

Returns:

The ratio of the applied change, of shape (...), equal to 1 where mask is False.

Return type:

torch.Tensor

property delay: int¶

Number of accepted updates not yet merged into the stored inverse.

Returns:

The number of delayed updates, between 0 and max_delay - 1.

Return type:

int

classmethod from_inverse(inverse: Tensor, *, pfaffian: Tensor | None = None, max_delay: int = 8) → PfaffianUpdater¶

Start from an inverse that is already known, without inverting a matrix or computing its Pfaffian.

Parameters:

• inverse – Inverse of the initial skew-symmetric matrix, of shape (..., n, n).

• pfaffian – Pfaffian of the initial matrix, of shape (...). When None, pfaffian starts at 1 and is the product of the accepted ratios, i.e. the Pfaffian relative to the start.

• max_delay – Number of accepted updates after which they are merged into the inverse, at least 1.

Returns:

The updater.

Return type:

PfaffianUpdater

Raises:

• ValueError – If max_delay is smaller than 1 or inverse is not square.

• TypeError – If the dtype of inverse is not supported.

property inverse: Tensor¶

Inverse of the current matrix, with the delayed updates applied (the state is not changed).

Returns:

The inverse, of shape (..., n, n).

Return type:

torch.Tensor

property max_delay: int¶

Number of accepted updates after which they are merged into the stored inverse.

Returns:

The maximum delay.

Return type:

int

merge() → None¶

Fold the delayed updates into the stored inverse, A_0^{-1} <- A_tau^{-1}.

property pfaffian: Tensor¶

that of the matrix given to reset() times the accepted ratios.

Returns:

The Pfaffian, of shape (...).

Return type:

torch.Tensor

Type:

Pfaffian of the current matrix

ratio(update: Tensor) → Tensor¶

Ratio Pf(A - u J u^T) / Pf(A) for the current matrix A, without changing the state.

Parameters:

update – Update vectors u, of shape (..., n, 2k).

Returns:

The ratio, of shape (...).

Return type:

torch.Tensor

ratio_replace(index: Tensor | int, change: Tensor) → Tensor¶

Ratio for adding change to row index (and subtracting it from the column), without changing the state. See ratio_replace().

Parameters:

• index – Index of the row and column, an integer or an integer tensor of shape (...).

• change – Change of the row, of shape (..., n).

Returns:

The ratio, of shape (...).

Return type:

torch.Tensor

reset(matrix: Tensor) → None¶

Restart from matrix: recompute its inverse and Pfaffian and drop the delayed updates.

Useful to bound the round-off accumulated over many updates, with the current matrix kept by the caller.

Parameters:

matrix – Invertible skew-symmetric matrix of shape (..., n, n).

Raises:

• ValueError – If matrix is not square.

• TypeError – If the dtype of matrix is not supported.

torch_pfaffian.low_rank.skew_gradient module¶

class torch_pfaffian.low_rank.skew_gradient.SkewGradient(*args, **kwargs)¶

Bases: Function

Identity whose backward pass projects the gradient onto skew-symmetric matrices, (g - g^T) / 2.

A function of a skew-symmetric matrix only depends on it through its skew-symmetric directions, but how an implementation reads the matrix (an upper triangle, a transposed product) decides where autograd puts the gradient of a raw leaf matrix. pfaffian() returns the skew-split convention Pf(A) A^{-T} / 2; wrapping the skew-symmetric inputs of the low-rank updates in this identity gives them the same convention. The gradient of any parameter the matrix is built from (keeping it skew-symmetric) is unchanged, since only the skew-symmetric part of the gradient enters it.

The forward pass returns a view of the input and costs nothing; the backward pass costs one pass over the n x n gradient. The input is a matrix of shape (..., n, n); complex matrices use the transpose, not the conjugate transpose.

static apply_if_needed(matrix: Tensor) → Tensor¶

Wrap matrix in the identity only when a gradient will flow through it, so that inference pays nothing.

Parameters:

matrix – Skew-symmetric matrix of shape (..., n, n).

Returns:

matrix, wrapped when it requires a gradient and gradients are enabled.

Return type:

torch.Tensor

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

Skew-symmetric part of the incoming gradient.

Parameters:

• ctx – Unused context.

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

Returns:

(grad_output - grad_output^T) / 2.

Return type:

torch.Tensor

static forward(matrix: Tensor) → Tensor¶

static setup_context(ctx: FunctionCtx, inputs: tuple, output: Tensor) → None¶

Module contents¶