{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5024a218",
   "metadata": {},
   "source": [
    "# The `pfaffian` function\n",
    "\n",
    "The **Pfaffian** of a $2n \\times 2n$ skew-symmetric matrix $A$ (a square matrix with\n",
    "$A^{\\top} = -A$) is a polynomial in its entries, with several equivalent definitions.\n",
    "\n",
    "**Formal definition (sum over perfect matchings).** Let $\\Pi$ be the set of all *perfect matchings*\n",
    "of $\\{1, 2, \\dots, 2n\\}$: partitions into $n$ unordered pairs\n",
    "$\\alpha = \\{(i_1, j_1), \\dots, (i_n, j_n)\\}$ with every $i_k < j_k$. Then\n",
    "\n",
    "$$\\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},$$\n",
    "\n",
    "where $\\operatorname{sgn}(\\alpha)$ is the sign of the permutation that sends $(1, 2, \\dots, 2n)$ to\n",
    "$(i_1, j_1, \\dots, i_n, j_n)$. Equivalently, as a sum over the full symmetric group $S_{2n}$,\n",
    "\n",
    "$$\\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)}.$$\n",
    "\n",
    "**Recursive definition (expansion along a row).** The Pfaffian also satisfies a Laplace-like\n",
    "recursion. Writing $A_{\\hat{1}\\hat{\\jmath}}$ for the $(2n-2) \\times (2n-2)$ matrix obtained by\n",
    "deleting rows and columns $1$ and $j$ from $A$,\n",
    "\n",
    "$$\\operatorname{Pf}(A) = \\sum_{j=2}^{2n} (-1)^{j}\\, a_{1 j}\\, \\operatorname{Pf}\\!\\big(A_{\\hat{1}\\hat{\\jmath}}\\big),$$\n",
    "\n",
    "with the base cases $\\operatorname{Pf} = 1$ for the empty $0 \\times 0$ matrix and\n",
    "$\\operatorname{Pf}\\!\\begin{pmatrix} 0 & a \\\\ -a & 0 \\end{pmatrix} = a$ for the $2 \\times 2$ case.\n",
    "\n",
    "**Relation to the determinant.** All of these agree, and the Pfaffian is the *signed square root* of\n",
    "the determinant:\n",
    "\n",
    "$$\\operatorname{Pf}(A)^2 = \\det(A).$$\n",
    "\n",
    "Where $\\sqrt{\\det A}$ loses the sign, the Pfaffian keeps it. `torch_pfaffian.pfaffian` computes it for\n",
    "PyTorch tensors. This tutorial shows the inputs the function accepts and the outputs it returns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4eb7785e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:35.292871Z",
     "iopub.status.busy": "2026-07-13T21:25:35.291863Z",
     "iopub.status.idle": "2026-07-13T21:25:38.342053Z",
     "shell.execute_reply": "2026-07-13T21:25:38.341040Z"
    }
   },
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "from torch_pfaffian import pfaffian\n",
    "\n",
    "# Seed the RNG so the random example matrices below are reproducible.\n",
    "_ = torch.manual_seed(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80e28cd8",
   "metadata": {},
   "source": [
    "## Input: a single skew-symmetric matrix\n",
    "\n",
    "The basic input is a skew-symmetric matrix of shape `(2n, 2n)`, and the output is a scalar tensor:\n",
    "its Pfaffian. For a $2 \\times 2$ matrix $\\begin{pmatrix} 0 & a \\\\ -a & 0 \\end{pmatrix}$ the Pfaffian\n",
    "is simply $a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "794e34ac",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:38.345122Z",
     "iopub.status.busy": "2026-07-13T21:25:38.345122Z",
     "iopub.status.idle": "2026-07-13T21:25:38.357785Z",
     "shell.execute_reply": "2026-07-13T21:25:38.356776Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input shape : (2, 2)\n",
      "output      : tensor(-3.)\n"
     ]
    }
   ],
   "source": [
    "matrix = torch.tensor([[0.0, -3.0], [3.0, 0.0]])\n",
    "pf = pfaffian(matrix)\n",
    "\n",
    "print(\"input shape :\", tuple(matrix.shape))\n",
    "print(\"output      :\", pf)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03059147",
   "metadata": {},
   "source": [
    "Any skew-symmetric matrix works. A quick way to build one is `A = M - M.T`, which is\n",
    "skew-symmetric by construction. The defining identity $\\operatorname{Pf}(A)^2 = \\det(A)$ then\n",
    "holds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "10ebf677",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:38.361811Z",
     "iopub.status.busy": "2026-07-13T21:25:38.360301Z",
     "iopub.status.idle": "2026-07-13T21:25:38.368938Z",
     "shell.execute_reply": "2026-07-13T21:25:38.368938Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pf(A)      = 0.682274341583252\n",
      "Pf(A) ** 2 = 0.46549826860427856\n",
      "det(A)     = 0.4654996395111084\n"
     ]
    }
   ],
   "source": [
    "m = torch.randn(6, 6)\n",
    "a = m - m.T  # skew-symmetric: a.T == -a\n",
    "\n",
    "pf = pfaffian(a)\n",
    "print(\"Pf(A)      =\", pf.item())\n",
    "print(\"Pf(A) ** 2 =\", (pf**2).item())\n",
    "print(\"det(A)     =\", torch.linalg.det(a).item())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a46a1565",
   "metadata": {},
   "source": [
    "## Input: a batch of matrices\n",
    "\n",
    "Leading dimensions are treated as a batch. An input of shape `(..., 2n, 2n)` returns one Pfaffian\n",
    "per matrix, with shape `(...,)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bcfda368",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:38.372977Z",
     "iopub.status.busy": "2026-07-13T21:25:38.372977Z",
     "iopub.status.idle": "2026-07-13T21:25:38.380618Z",
     "shell.execute_reply": "2026-07-13T21:25:38.380113Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input shape : (4, 8, 8)\n",
      "output shape: (4,)\n",
      "output      : tensor([-22.6322,  11.7543,   5.5873, -16.7749])\n"
     ]
    }
   ],
   "source": [
    "batch = torch.randn(4, 8, 8)\n",
    "batch = batch - batch.transpose(-1, -2)  # skew-symmetric in the last two dims\n",
    "\n",
    "pfs = pfaffian(batch)\n",
    "print(\"input shape :\", tuple(batch.shape))\n",
    "print(\"output shape:\", tuple(pfs.shape))\n",
    "print(\"output      :\", pfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bdd646c",
   "metadata": {},
   "source": [
    "## Input dtype: real or complex\n",
    "\n",
    "The Pfaffian is defined for complex skew-symmetric matrices too. The output matches the input's\n",
    "dtype: a real input gives a real Pfaffian, a complex input gives a complex one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7f230a25",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:38.384451Z",
     "iopub.status.busy": "2026-07-13T21:25:38.383930Z",
     "iopub.status.idle": "2026-07-13T21:25:38.392680Z",
     "shell.execute_reply": "2026-07-13T21:25:38.391647Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "real input    -> torch.float32 | -3.7162179946899414\n",
      "complex input -> torch.complex64 | (-0.5062580108642578+0.7343237996101379j)\n"
     ]
    }
   ],
   "source": [
    "real = torch.randn(4, 4)\n",
    "real = real - real.T\n",
    "\n",
    "complex_matrix = torch.randn(4, 4, dtype=torch.complex64)\n",
    "complex_matrix = complex_matrix - complex_matrix.transpose(-1, -2)\n",
    "\n",
    "print(\"real input    ->\", pfaffian(real).dtype, \"|\", pfaffian(real).item())\n",
    "print(\"complex input ->\", pfaffian(complex_matrix).dtype, \"|\", pfaffian(complex_matrix).item())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f1d6ffe",
   "metadata": {},
   "source": [
    "## Output: signed value or magnitude\n",
    "\n",
    "By default `pfaffian` returns the **signed** Pfaffian. Pass `sign=False` to get only its magnitude\n",
    "$|\\operatorname{Pf}(A)|$, which takes a cheaper determinant-based path when the sign is not needed.\n",
    "Here the $2 \\times 2$ matrix has a negative Pfaffian, so the two differ by the sign:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0ec435e8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:38.395300Z",
     "iopub.status.busy": "2026-07-13T21:25:38.395300Z",
     "iopub.status.idle": "2026-07-13T21:25:38.401040Z",
     "shell.execute_reply": "2026-07-13T21:25:38.401040Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "signed (default) : -3.0\n",
      "magnitude        : 3.0\n"
     ]
    }
   ],
   "source": [
    "print(\"signed (default) :\", pfaffian(matrix).item())\n",
    "print(\"magnitude        :\", pfaffian(matrix, sign=False).item())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4fd3feb",
   "metadata": {},
   "source": [
    "## Output: differentiable\n",
    "\n",
    "If the input requires gradients, so does the output: `pfaffian` is fully compatible with\n",
    "`torch.autograd`, so it composes with any PyTorch computation and backpropagates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "08bba05a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-13T21:25:38.403668Z",
     "iopub.status.busy": "2026-07-13T21:25:38.403668Z",
     "iopub.status.idle": "2026-07-13T21:25:38.409793Z",
     "shell.execute_reply": "2026-07-13T21:25:38.409266Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input.requires_grad : True\n",
      "gradient shape      : (6, 6)\n"
     ]
    }
   ],
   "source": [
    "a_grad = a.clone().requires_grad_(True)\n",
    "pfaffian(a_grad).backward()\n",
    "\n",
    "print(\"input.requires_grad :\", a_grad.requires_grad)\n",
    "print(\"gradient shape      :\", tuple(a_grad.grad.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69e60b67",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "| Input to `pfaffian` | Output |\n",
    "| --- | --- |\n",
    "| skew-symmetric matrix `(2n, 2n)` | scalar Pfaffian |\n",
    "| batch `(..., 2n, 2n)` | Pfaffians of shape `(...,)` |\n",
    "| real or complex dtype | Pfaffian of the same dtype |\n",
    "| `requires_grad=True` | differentiable output |\n",
    "| CPU or CUDA tensor | result on the same device |\n",
    "\n",
    "In every case the output shares the input's dtype, device, and backend, and returns a value whose\n",
    "square is $\\det(A)$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}