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.

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
• What problem this solves
• New in 1.2: blocks, adjoints and expressions
• Tutorial: a prepare-and-measure scenario
• Building the SDP
• Localizing matrices and sum rules
• Block moment matrices
• API reference
• Performance
• Upgrading from 0.x

---

Installation

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

From a checkout:

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:

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:

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, 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:

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; the full list is in 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.

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)

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:

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:

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.

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

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.

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:

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

4. Build

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 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:

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

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

Index the model by monomial or by variable index:

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:

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:

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.

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:

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

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:

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

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

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:

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).

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).

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.

---

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.

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 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.