Metadata-Version: 2.4
Name: MoMPy
Version: 1.2.0
Summary: Moment matrices for SDP hierarchy relaxations
Author-email: Carles Roch i Carceller <chalswater@gmail.com>
License: MIT
Project-URL: Homepage, https://github.com/chalswater/MoMPy
Project-URL: Repository, https://github.com/chalswater/MoMPy
Project-URL: Issues, https://github.com/chalswater/MoMPy/issues
Keywords: semidefinite programming,SDP,moment matrix,NPA hierarchy,quantum information,noncommutative polynomial optimization
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Physics
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.20
Provides-Extra: cvxpy
Requires-Dist: cvxpy>=1.3; extra == "cvxpy"
Provides-Extra: test
Requires-Dist: pytest>=7; extra == "test"
Requires-Dist: cvxpy>=1.3; extra == "test"
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Requires-Dist: cvxpy>=1.3; extra == "dev"
Requires-Dist: ruff; extra == "dev"
Requires-Dist: build; extra == "dev"
Dynamic: license-file

# MoMPy

**Moment matrices for SDP hierarchy relaxations.**

MoMPy builds the moment matrix of a semidefinite relaxation and works out, for
you, which of its entries are forced to be equal or zero by the algebraic
properties of your operators — rank-1 projectors, orthogonal measurements,
commutation. You describe the operators; MoMPy hands back a matrix of SDP
variable indices ready to drop into CVXPY. Since 1.2 operators may also be
non-Hermitian, obey arbitrary substitution rules such as the matrix units of a
subsystem, and be combined into polynomials that the CVXPY helper turns into
moments, localizing matrices and relations.

```python
from MoMPy import OperatorSet, MomentProblem

ops = OperatorSet()
R = ops.add_family(3, idempotent=True)     # three pure states
M = ops.add_povm_family(2, 2)              # M[y][b]: two binary measurements
ops.declare_commuting(R, R)                # the states commute with each other

monomials  = list(R) + [m for row in M for m in row]
monomials += [[R[x], M[y][b]] for x in range(3) for y in range(2) for b in range(2)]

mm = MomentProblem(monomials, ops.algebra(), dim=1).build()
print(mm.summary())
```

```
MomentMatrix: 20 x 20 (19 monomials + identity)
  block size (dim)   : 1
  SDP variables      : 64
  compression        : 400 entries -> 64 variables (6.2x)
  zero entries       : 64
  distinct words seen: 885
  build time         : 0.012 s
```


> ### `cyclicity`: tracial or state moments? Read this before your first build.
>
> `MomentProblem(..., cyclicity=True)` (the default) uses **tracial** moments
> `Tr(u v†)`, which are cyclic. `cyclicity=False` uses **state** moments
> `<psi|u v†|psi>`, which are not.
>
> Cyclicity is valid when your figure of merit really is a trace with the state
> *inside* the algebra — prepare-and-measure scenarios, `Tr(rho_x M_b)`. It is
> **not** valid for Bell/NPA problems. Imposing it there over-constrains the
> program: CHSH at level 1+AB returns 2.0000 instead of Tsirelson's 2.8284, so
> it is not an upper bound on the quantum value at all. At level 1 both agree,
> which makes the error easy to miss.
>
> | Problem | Use |
> |---|---|
> | Bell, NPA, device-independent | `cyclicity=False` |
> | Prepare-and-measure, dimension witnesses | `cyclicity=True` (the default) |
> | Unsure | `cyclicity=False` (fewer relations, so never invalid) |

---

## Contents

- [Installation](#installation)
- [What problem this solves](#what-problem-this-solves)
- [New in 1.2: blocks, adjoints and expressions](#new-in-12-blocks-adjoints-and-expressions)
- [Tutorial: a prepare-and-measure scenario](#tutorial-a-prepare-and-measure-scenario)
- [Building the SDP](#building-the-sdp)
- [Localizing matrices and sum rules](#localizing-matrices-and-sum-rules)
- [Block moment matrices](#block-moment-matrices)
- [API reference](#api-reference)
- [Performance](#performance)
- [Upgrading from 0.x](#upgrading-from-1x)

---

## Installation

```bash
pip install MoMPy            # core, needs only numpy
pip install MoMPy[cvxpy]     # plus the CVXPY helpers
```

From a checkout:

```bash
pip install -e ".[dev]"
pytest
```

---

## What problem this solves

Take a prepare-and-measure scenario. Alice encodes a message `x` in a quantum
state `R[x]` and sends it to Bob, who measures with `M[y][b]` and observes `b`.
The observable statistics are `p(b|x,y) = Tr(R[x] @ M[y][b])`, and you want to
maximise some linear functional of them over *all* states and measurements.

That optimisation is not an SDP. The standard relaxation makes it one: list
monomials in your operators, `L = {1, R[x], M[y][b], R[x] R[x'], R[x] M[y][b], ...}`,
and form the matrix `G[u,v] = Tr(u v†)` over `u, v ∈ L`. `G` is positive
semidefinite by construction and your objective lives inside it, so maximising
over PSD `G` gives an upper bound.

The tedious part is that many entries of `G` are secretly the same variable.
If `R[x]` is a pure state then `Tr(R[x])` and `Tr(R[x] R[x])` are equal. If
`M[y][b]` is a projective measurement then `Tr(R[x] M[y][0] M[y][1])` is
identically zero. Miss these identifications and your relaxation is looser than
it should be; get them wrong and it is not a valid bound at all.

MoMPy finds them. You declare the properties, it computes the equivalence
classes and returns the matrix.

**Applicable to** any optimisation expressible as an SDP relaxation over traces
of operator monomials: NPA / device-independent bounds, prepare-and-measure
scenarios, dimension witnesses, randomness certification, joint measurability.

---

## New in 1.2: what you can now describe

Up to 1.1, MoMPy could express three kinds of relation — idempotency,
orthogonality, commutation — and assumed, everywhere and silently, that every
operator is Hermitian. That covers Bell scenarios and many prepare-and-measure
ones, and nothing about it changes.

1.2 widens what a hierarchy can be *about*. Four additions, each independent
of the others:

**Operators need not be Hermitian.** Declare the involution and every part of
the engine uses the true adjoint `(l₁…l_n)† = l_n†…l₁†` instead of plain
reversal:

```python
ops.declare_adjoint(a, b)                        # a† = b
A = ops.add_tensor(k, k, adjoint="transpose")    # A_ij† = A_ji
```

This is a correctness matter, not a convenience: the `hermitian` move
identifies a word with its reversal, which is the adjoint *only* if every
letter is Hermitian. Put a non-Hermitian letter in the alphabet without
declaring it and the relaxation silently stops being a valid bound. Useful
whenever the alphabet contains isometries, unitaries, Kraus operators,
annihilation/ladder operators, matrix units, or the off-diagonal blocks of a
Hermitian operator.

**Relations need not be idempotency or orthogonality.** Any word can be
declared to reduce to another word, to the identity, or to zero:

```python
ops.declare_substitution((U_dag, U), ())        # an isometry: U†U = 1
ops.declare_substitution((V, V, V), ())         # a cube root: V³ = 1
ops.declare_substitution((E01, E10), E00)       # matrix units
ops.declare_substitution((P, F, Q), 0)          # a vanishing sandwich
```

Patterns of any length, matched at every position. This is what lets group
algebras, partial isometries, mutually unbiased structures, Naimark blocks and
subsystem decompositions be described directly rather than approximated by
extra constraints.

**Operator inequalities and sum rules get first-class support.** `0 ≤ M ≤ 1`,
`ρ ≥ 0`, `Σ_b Π_b = 1`, `N² = N` are not relations between *words*, so they
cannot be declarations at all. The CVXPY model now builds them — localizing
matrices for the inequalities, contextual moment equalities for the sums — see
[Localizing matrices and sum rules](#localizing-matrices-and-sum-rules), which
is probably the single biggest addition in this release.

**Operators can be written as operators.** Labels multiply and add into
polynomials, so the object you are constraining is the object you type:

```python
N = sum(E[i][j] * A[i][j] for i in range(k) for j in range(k)) + F
ct += model.localizing(ONE - N, mons, weight=R[x])
ct += model.relation(N * N - N, contexts=2, weight=R[x])
```

Three smaller items: `state_monomials` for prepare-and-measure hierarchies
that keep the states in the algebra, `dedupe="operators"` to drop monomials
that are redundant as operators, and lookups that reduce a word before
reporting it missing.

### Which tier does a given relation belong to?

| Your relation | Where it goes |
|---|---|
| One word equals another word, or zero (`P P = P`, `U†U = 1`, `E_ij E_jm = E_im`) | A **declaration** — it merges SDP variables, so it costs nothing and tightens everything |
| A sum or a scalar multiple (`Σ_b Π_b = 1`, `N² = N`, `A P A = ½ A`) | `model.relation(...)` — moment equalities, imposed in every context |
| An inequality / positivity (`0 ≤ M ≤ 1`, `ρ ≥ 0`, `X ≥ 0`) | `model.localizing(...)` |

The first tier is strictly the strongest, so push a relation as far up this
table as it will go. A declaration cannot carry a coefficient or a sum — a
pattern maps to *one* word or to zero — which is exactly where the second tier
takes over.

A complete worked problem using all of it is in
[`examples/energy_constrained_discrimination_example.py`](examples/energy_constrained_discrimination_example.py);
the full list is in [`PATCH_NOTES_1.2.0.md`](PATCH_NOTES_1.2.0.md).

---

## Tutorial: a prepare-and-measure scenario

### 1. Allocate operators

Operators are integer labels. `OperatorSet` allocates them and remembers their
properties, so you never keep a counter by hand. **Label `0` is reserved for
the identity** and is added to the matrix automatically.

```python
from MoMPy import OperatorSet

nX, nY, nB = 3, 2, 2

ops = OperatorSet()
R = ops.add_family(nX, idempotent=True)     # R[x],   pure states
M = ops.add_povm_family(nY, nB)             # M[y][b], projective measurements
```

`add_povm_family` registers each measurement's outcomes as an orthogonal set
*and* as projectors, which is the usual projective assumption. Override with
`add_povm(n, idempotent=False, orthogonal=False)` if you need something else.

### 2. Declare the relations

Relations declared here are folded into the moment matrix itself: words related
by them share one SDP variable, which shrinks the program and tightens the
relaxation at no cost.

| Relation | Meaning | How to declare |
|---|---|---|
| **Idempotent** | `P @ P == P` | `add_family(..., idempotent=True)` or `ops.declare_idempotent([...])` |
| **Orthogonal** | `P_i @ P_j == 0` for `i != j` | `add_povm(...)` or `ops.declare_orthogonal([...])` |
| **Orthogonal, two groups** *(1.2)* | `a @ b == 0` for every `a in A`, `b in B` | `ops.declare_orthogonal(A, B)` |
| **Commuting** | `a @ b == b @ a` | `ops.declare_commuting(A, B)` |
| **Adjoint** *(1.2)* | `a† == b`; anything undeclared is Hermitian | `ops.declare_adjoint(a, b)`, or `add_tensor(..., adjoint="transpose")` |
| **Substitution** *(1.2)* | a word reduces to another word, to `1`, or to `0` | `ops.declare_substitution(pattern, replacement)` |
| **Matrix units** *(1.2)* | the whole table `E_ij E_lm == δ_jl E_im`, plus adjoints | `ops.declare_matrix_units(E)` |

```python
ops.declare_commuting(R, R)                 # every R[x] commutes with every R[x']
```

`declare_commuting(A, B)` means *every* label in `A` commutes with *every* label
in `B`. Pass the same list twice for "all of these commute with each other".
`declare_orthogonal` follows the same convention: one group means "mutually
orthogonal", two groups mean "each of these against each of those".

**Adjoints.** Every label is Hermitian unless you say otherwise, so scenarios
written before 1.2 need no change. When an operator is *not* Hermitian — an
isometry, a Kraus operator, a matrix unit `|i><j|`, an off-diagonal block —
declare its partner:

```python
ops.declare_adjoint(U, U_dag)                   # U† = U_dag

A = ops.add_tensor(k, k, adjoint="transpose")   # A[i][j]† = A[j][i]
```

`adjoint="transpose"` is the shorthand for a square array whose entries behave
like the entries of a matrix under the dagger: the adjoint of `A[i][j]` is
`A[j][i]`, and the diagonal is Hermitian. That is the structure of any
Hermitian operator written in blocks, `A = Σ_ij |i><j| ⊗ A_ij`. Leading
dimensions are untouched, so a stack of block matrices is
`add_tensor(n_b, k, k, adjoint="transpose")`. For anything more exotic, pass a
callable mapping an index tuple to the index tuple of its adjoint.

**Substitutions.** A pattern of any length, and a right-hand side that is a
word, `()` for the identity, or `0` for zero:

```python
ops.declare_substitution((U_dag, U), ())        # U†U = 1: an isometry
ops.declare_substitution((V, V, V), ())         # V³ = 1
ops.declare_substitution((a, b), c)             # a b = c
ops.declare_substitution((P, F, Q), 0)          # P F Q = 0
```

A pattern reduces to a single word or to zero — it cannot produce a sum or
carry a coefficient. Relations like `Σ_b Π_b = 1` or `A P A = ½A` are moment
relations instead; see
[Localizing matrices and sum rules](#localizing-matrices-and-sum-rules).

**Matrix units.** The relations of a subsystem's basis, `E_ij = |i><j|`, come
as a set:

```python
E = ops.add_tensor(k, k)
ops.declare_matrix_units(E)      # E_ij E_lm = δ_jl E_im, E_ij† = E_ji, E_ii² = E_ii
```

This is the standard way to expose one subsystem's degrees of freedom inside a
hierarchy: resolve an operator that acts jointly on two systems in a basis of
the first, and the blocks acting on the second appear explicitly and commute
with the units.

### 3. Choose your monomials

The hierarchy level is just which monomials you include. Longer words give a
tighter bound and a bigger matrix.

```python
monomials  = list(R)                                      # first order
monomials += [m for row in M for m in row]
monomials += [[R[x], M[y][b]]                             # second order
              for x in range(nX) for y in range(nY) for b in range(nB)]
monomials += [[R[x], R[xx], R[xxx]]                       # some third order
              for x in range(nX) for xx in range(nX) for xxx in range(nX)]
```

A monomial is a bare label or a list of labels read left to right as a product.
For the standard "all words up to length k" there is a shortcut:

```python
from MoMPy import generate_monomials
monomials = generate_monomials(list(R) + flat_M, level=2)
```

### 4. Build

```python
from MoMPy import MomentProblem

mm = MomentProblem(monomials, ops.algebra(), dim=1).build(progress=True)
```

`dim` is the one parameter with no default: it is the side length of the
block that will back each entry once you reach `to_cvxpy` (see
[Block moment matrices](#block-moment-matrices) below). `dim=1` is the
ordinary scalar moment matrix used throughout this tutorial section.

`mm.matrix` is an integer NumPy array: `mm.matrix[r, c]` is the index of the SDP
variable at that position. Equal indices mean the same variable.

Look up the variable for any monomial:

```python
mm.index_of([R[0], M[1][0]])     # the variable holding Tr(R0 M10)
mm.identity_index                # the variable holding Tr(1)
mm.zero_index                    # the class of monomials forced to zero
mm.equivalents([R[0]])           # every monomial equal to Tr(R0)
```

---

## Building the SDP

### With the CVXPY helper

```python
model = mm.to_cvxpy()
ct = list(model.constraints)          # G >> 0, and zeros pinned to zero
```

Index the model by monomial or by variable index:

```python
model[[R[0], M[1][0]]]     # scalar expression for Tr(R0 M10)
model.identity             # Tr(1)
```

> **`Tr(1)` is the dimension, not 1.** In a tracial relaxation the identity
> variable equals the Hilbert-space dimension. MoMPy deliberately does *not*
> constrain it. Add `ct.append(model.identity == 1)` only if you are using the
> state-vector NPA convention where moments are `<psi| w |psi>`.

### Normalisation constraints

`sum_b M[y][b] == 1` is a linear relation between variables, so it must be added
to the program. MoMPy finds every place it applies:

```python
for y in range(nY):
    ct += model.apply(mm.normalisation_constraints(M[y]))
```

For joint measurability, where a parent POVM marginalises onto a single
operator:

```python
ct += model.apply(mm.marginal_constraints(joint=B_labels, marginal=M[0][0]))
```

### Localizing matrices and sum rules

*New in 1.2, and the largest addition in the release.*

A moment matrix can only talk about *moments*. Two very common kinds of
statement are not moments, and before 1.2 you had to build both by hand:

- **Inequalities.** `0 ≤ M ≤ 1`, `ρ ≥ 0`, `1 - N ≥ 0`. These are operator
  positivity statements, not numbers.
- **Sum rules.** `Σ_b Π_b = 1`, `N² = N`, `A_ij = Σ_l A_il A_lj`. These involve
  a sum, so they cannot be word rewrites either.

#### Inequalities: `model.localizing`

For a positive operator `X ≥ 0` and any list of monomials `u`, the matrix

```
L[u, v] = Tr( u† X v )
```

is positive semidefinite, because it is a Gram matrix of the vectors
`√X u`. Constraining `L >> 0` is what makes the rest of the program aware that
`X` is positive; it is called a *localizing matrix*, and it is the standard
device for inequality constraints in noncommutative polynomial optimisation.

```python
mons = [ONE] + letters                      # the longer the list, the tighter

ct += model.localizing(N, mons)             # N >= 0
ct += model.localizing(ONE - N, mons)       # N <= 1
```

`X` is any polynomial, so `ONE - N` and `N * N` work as readily as a bare
label. The monomial list is yours to choose: `[ONE]` alone gives just
`Tr(X) ≥ 0`, and each monomial you add strengthens the constraint at the cost
of a larger block. The words `u† X v` must lie within the hierarchy, so the
useful rule of thumb is `deg(u) + deg(X) + deg(v) ≤ 2 × level`.

**The `weight` argument.** Passing an operator multiplies from the left:

```python
ct += model.localizing(ONE - N, mons, weight=R[x])     # Tr(ρ_x u† (1-N) v) >> 0
```

This is a *state* localizing matrix, and in a tracial hierarchy — where the
states live in the algebra and `Tr(1)` is a free variable standing for the
Hilbert-space dimension — it is usually the version that carries the physics.
The unweighted moments `Tr(u† X v)` are inner products in an unbounded space
and constrain very little on their own; weighting by a state is what ties `X`
to the actual preparation. If a tracial relaxation comes back trivially loose,
a missing weight is the first thing to check.

Use `model.localizing_matrix(...)` if you want the expression itself, to
inspect it or to read its value after solving.

#### Sum rules: `model.relation`

```python
ct += model.relation(N * N - N, contexts=2)                    # N² = N
ct += model.relation(sum(Pi) - ONE, contexts=1)                # Σ_b Π_b = 1
ct += model.relation(A[i][j] - sum(A[i][l] * A[l][j] for l in range(k)),
                     contexts=[ONE] + letters, weight=R[x])
```

Imposing `Tr(poly) == 0` alone is far weaker than imposing `Tr(u poly v) == 0`
for every context the matrix can express, so `relation` does the latter.
`contexts` is a list of monomials, or an integer `d` meaning every monomial of
the hierarchy of length at most `d` plus the identity. Contexts whose words
fall outside the hierarchy are skipped — dropping constraints only relaxes the
program, so the bound stays valid — and `strict=True` raises instead.

#### Both return lists

`localizing` and `relation` both return **lists of constraints**, so they read
the same way at the call site:

```python
ct = list(model.constraints)                     # structural: PSD, zeros
ct += model.localizing(ONE - N, mons, weight=R[0])
ct += model.relation(N * N - N, contexts=2, weight=R[0])
ct += model.apply(mm.normalisation_constraints(M[0]))
ct += [model[R[0] * E[0][0]] >= 1 - omega]       # a plain CVXPY comparison
```

The last line is bracketed because it is an ordinary CVXPY constraint, as any
hand-written comparison is; everything MoMPy generates comes back as a list.

### Problem-specific constraints

```python
ct += [model[[R[x]]] == 1.0 for x in range(nX)]              # states are normalised
ct += [model[[R[x], R[xx]]] >= d for x in range(nX) for xx in range(nX)]
```

### Solve

```python
import cvxpy as cp

W = sum(model[[R[x], M[0][x]]] for x in range(nX))
problem = cp.Problem(cp.Maximize(W), ct)
problem.solve(solver=cp.SCS)
print(problem.value)
```

Any SDP solver works — SCS and Clarabel ship with CVXPY; MOSEK is free with an
academic licence.

### Without CVXPY

Nothing ties you to CVXPY. Allocate one variable per index and read the matrix:

```python
variables = {i: make_variable() for i in mm.variable_indices}
variables[mm.zero_index] = 0.0
G = [[variables[mm.matrix[r, c]] for c in range(mm.n)] for r in range(mm.n)]
```

Constraint objects expose plain integers via `.lhs` and `.rhs`, so
`mm.normalisation_constraints(...)` is usable with any modelling layer.

---

## Block moment matrices

Set `dim=d` for `d > 1` and every entry of the matrix becomes a `d x d` block
instead of a scalar — for relaxations whose "moments" are themselves
operators on a `d`-dimensional Hilbert space, rather than numbers. Cyclicity
practically never holds for these: `u v` and `v u` are genuinely different
blocks, so `cyclicity=False` is the right choice for essentially every block
hierarchy, and `hermitian=False` too whenever a block and its adjoint are
meant to be different blocks (the general case).

```python
bm = MomentProblem(monomials, ops.algebra(), dim=d, cyclicity=False, hermitian=False).build()
model = bm.to_cvxpy()            # dim x dim CVXPY blocks, read straight off bm.dim
model[[R[0], M[1][0]]]           # a dim x dim expression, not a scalar
```

`to_cvxpy` is the same function used for scalar matrices above — it reads
`matrix.dim` and builds scalars or blocks accordingly, so there is nothing
extra to call or import for the block case. Everything else (constraints,
`.apply()`, `.normalisation_constraints()`, indexing by monomial or by
variable index) works exactly as in the scalar walkthrough above.

---

## API reference

### Describing a problem

| Object | Purpose |
|---|---|
| `OperatorSet` | Allocates labels, records properties, emits an `Algebra` |
| `Algebra(idempotents, orthogonal_sets, commuting_pairs, adjoint, substitutions)` | The relations, if you prefer to build them by hand |
| `generate_monomials(letters, level)` | All words up to a given length |
| `state_monomials(states, letters, level)` | The words `ρ_x w` (1.2) |
| `Poly`, `Label`, `ONE`, `op` | Operator expressions (1.2) |
| `MomentProblem(monomials, algebra, *, dim, cyclicity=True, hermitian=True, dedupe=True)` | One class for every relaxation: scalar or block, tracial or state |
| `MomentProblem.from_levels(letters, level, extra=..., states=..., state_level=..., dim=...)` | Shortcut constructor |

`OperatorSet` declarations: `declare_idempotent`, `declare_orthogonal` (one
group, or two groups since 1.2), `declare_commuting`, and — new in 1.2 —
`declare_adjoint(a, b)`, `declare_substitution(pattern, replacement)`,
`declare_matrix_units(E)`.

One class covers what used to be four: `MomentProblem(m, a, dim=1)` is the
tracial relaxation `Tr(u v†)`; add `cyclicity=False` for state moments
`<psi|u v†|psi>` — **use this for NPA/Bell** — and `dim=d>1` for a block
hierarchy whose entries are `d x d` operators (see
[Block moment matrices](#block-moment-matrices)).

### `MomentProblem.build(progress=False)` → `MomentMatrix`

| Attribute | Meaning |
|---|---|
| `.matrix` | `(n, n)` integer array of variable indices |
| `.n`, `.shape` | Matrix size |
| `.monomials` | Generating monomials, excluding the identity |
| `.word_at(r, c)` | Explicit operator word behind an entry |
| `.words` | Full nested list of words (built lazily) |
| `.map_table` | `MapTable`: monomial → index |
| `.variable_indices`, `.n_variables` | The distinct variables present |
| `.zero_index`, `.identity_index` | Reserved classes |
| `.has_zeros` | Whether orthogonality forced anything to zero |
| `.stats` | Build diagnostics |
| `.index_of(w)`, `.get(w, default)` | Lookup; reduces the word first (1.2), `index_of` raises `UnknownMonomial` |
| `.state_block(state)` | The sub-block indexed by `ρ_x w` (1.2) |
| `.equivalents(w)` | All monomials sharing `w`'s variable |
| `.summary()` | Human-readable report |
| `.normalisation_constraints(povm)` | `sum(povm) == 1` constraints |
| `.marginal_constraints(joint, marginal)` | `sum(joint) == marginal` constraints |
| `.to_cvxpy(dim=None, psd=True, complex=None, normalise_identity=False)` | CVXPY model, scalar or block per `.dim` |
| `.to_legacy()` | The 0.x five-tuple |
| `.dim`, `.cyclicity`, `.hermitian` | The three flags the matrix was built with |

### Options

- **`dim`** — side length of the block each SDP variable becomes in
  `to_cvxpy`. No default: declare it explicitly, even as `dim=1` for an
  ordinary scalar matrix.
- **`cyclicity`** — identify each word with its cyclic rotations, i.e. treat
  an entry as a trace `Tr(u v)` rather than an operator product `u v`. Default
  `True`. See the callout above — `False` is what NPA/Bell problems need.
- **`hermitian`** — identify each word with its **adjoint**. For Hermitian
  operators this says the moment matrix is real symmetric, i.e. the variables
  are `Re Tr(w)`. Default `True`. Set `False` to build a complex Hermitian
  SDP, or for a block hierarchy where a block and its adjoint should be
  independent. Since 1.2 the adjoint respects `declare_adjoint`; with no
  declaration it is plain reversal, exactly as before.
- **`dedupe`** — drop repeated monomials, which only add linearly dependent
  rows and columns. Default `True`. Since 1.2, `"operators"` also drops
  monomials that vanish or that equal an earlier one as operators.

### The CVXPY model (1.2)

| Method | Purpose |
|---|---|
| `model[poly]`, `model.value(poly, strict=)` | Evaluate an operator polynomial |
| `model.localizing(poly, monomials, weight=)` | **List of constraints** making `poly >= 0` |
| `model.localizing_matrix(poly, monomials, weight=)` | The localizing matrix expression itself |
| `model.relation(poly, contexts=, weight=)` | **List of constraints** `Tr(u poly v) == 0` |

Both `localizing` and `relation` return lists, so `ct += model.localizing(...)`
and `ct += model.relation(...)` read identically. See
[Localizing matrices and sum rules](#localizing-matrices-and-sum-rules).

---

## Performance

Version 2 replaces the per-monomial linear scans with canonical tuple words, a
breadth-first closure that memoises every word it has already seen, and a
union-find over classes. Each distinct word is expanded exactly once for the
whole build, and monomial lookup is a dict probe rather than a scan over every
word in every class.

Measured on the scenarios in `examples/`:

| Scenario | Matrix | 0.x | 1.x | Speedup |
|---|---|---|---|---|
| NPA CHSH level 1 | 9×9 | 0.01 s | 0.007 s | ~1× |
| NPA CHSH level 1+AB | 25×25 | 0.02 s | 0.012 s | 2× |
| PAM dimension, 3rd order | 84×84 | 41.7 s | 0.041 s | **1027×** |
| PAM dimension, 2nd+3rd order | 105×105 | 52.4 s | 0.057 s | **919×** |
| PAM dimension, nX=4 | 137×137 | 528 s | 0.178 s | **2960×** |

Scaling is now roughly linear in the number of matrix entries:

| Scenario | Matrix | Entries | Variables | Time |
|---|---|---|---|---|
| NPA 3 settings, 3 outcomes, 1+AB | 100×100 | 10 000 | 1 370 | 0.47 s |
| NPA 5 settings, 3 outcomes, 1+AB | 256×256 | 65 536 | 11 237 | 3.7 s |
| PAM 6 states, order 3 | 287×287 | 82 369 | 381 | 0.86 s |
| PAM 8 states, order 3 | 639×639 | 408 321 | 1 670 | 5.8 s |

---

## Correctness

The equivalence classes are checked against a deliberately naive brute-force
closure oracle over 720 randomised scenarios, covering tracial and block modes
with and without reversal symmetry. The induced partitions match exactly.

On top of that, `tests/test_physics.py` solves real SDPs (CHSH → 2√2, a fully
commutative algebra → the local bound 2, state discrimination → 1) and plugs
explicit matrices in for the operator labels to confirm numerically that every
monomial sharing a variable really does have the same trace and that the zero
class really vanishes.

```bash
pytest                     # everything
pytest tests/test_api.py   # fast unit tests only
```

---

## Upgrading from 0.x

**Your existing scripts keep working.** `from MoMPy.MoM import *` still gives
you `MomentMatrix`, `fmap`, `normalisation_contraints` and friends, returning
the same five outputs.

Two fixes do change the numbers you get, both in the direction of a tighter and
more correct relaxation. See [`MIGRATION.md`](MIGRATION.md) for the details and
for how to port to the new API.

---

## Citing and contact

Author: Carles Roch i Carceller — <chalswater@gmail.com>
Repository: <https://github.com/chalswater/MoMPy> · MIT licence.
