PML-361: phase errors on Unitary-based circuit using a decomposition

[CassNot]

Summary

Currently, if the circuit we give the QuantumLayer is Unitary-based and we have a NoiseModel, the phase error cannot be applied.
This is problematic because:

• Real photonic hardware has phase imprecision everywhere (not just named parameters)
• Black-box unitaries represent physical device components that are subject to the same fabrication tolerances
• Simulations without this noise don't match real hardware behavior (like the ReservoirClassifier)

Here, I am adding automatic decomposition of black-box Unitary components into Clements MZI meshes when circuit phase noise is configured. This enables applying phase noise (phase_imprecision and phase_error) to every physical phase in the circuit, not just parameterized phase shifters.

We should now be able to model realistic hardware imperfections affecting all circuit phases, including those inside black-box unitary blocks (like ReservoirClassifier)

Related Issue

PML-361

Type of change

• Bug fix
• New feature
• Documentation update
• Refactor / Cleanup
• Performance improvement
• CI / Build / Tooling
• Breaking change (requires migration notes)

Proposed changes

New helper function _decompose_unitaries() that:

1. Detects black-box Unitaries: When phase noise is configured, the circuit is scanned for non-PERM Unitary components
2. Decomposes to MZI mesh: Each unitary is fitted to a rectangular Clements MZI mesh using Perceval's CircuitOptimizer
3. Preserves relative phase: The global phase difference is folded into the output phase-shifter layer
4. Enables noise application: The decomposed mesh structure makes every phase explicit and subject to configured noise

Circuit with black-box Unitary
            ↓
   (phase noise configured?)
    ↓              ↓
   YES            NO
    ↓              ↓
Decompose      Fast Path
  to MZI    (unchanged)
    ↓
MZI Mesh (fixed + phase shifters)
    ↓
Apply phase_imprecision & phase_error
    ↓
Noisy Unitary

How to test / How to run

1. Command lines

pytest tests/pcvl_pytorch/test_unitary_decomposition_noise.py

Performance considerations (optional)

Documentation

• User docs updated (Sphinx)
• Examples / notebooks updated
• Docstrings updated
• Updated the API

Checklist

• PR title includes Jira issue key (e.g., PML-126)
• "Related Jira ticket" section includes the Jira issue key (no URL)
• Code formatted (ruff format)
• Lint passes (ruff)
• Static typing passes (mypy) if applicable
• Unit tests added/updated (pytest)
• Tests pass locally (pytest)
• Tests pass on GPU (pytest)
• Test coverage not decreased significantly
• Docs build locally if affected (sphinx)
• With this command: SPHINXOPTS="-W --keep-going -n" make -C docs clean html the docs are built without any warning or errors.
• New public classes/methods/packages are added in the API following the methodology presented in other files.
• Dependencies updated (if needed) and pinned appropriately
• PR description explains what changed and how to validate it

[CassNot]

@claude review

[claude]

Claude finished @CassNot's task in 6m 54s —— View job

---

Code review — PML-361 (phase errors on Unitary-based circuits)

• Gather context (AGENTS.md, PR diff vs origin/release/0.4.1)
• Review _decompose_unitaries() and its integration in CircuitConverter
• Review new tests (test_unitary_decomposition_noise.py, test_noise_contract.py, test_reservoir_classifier.py)
• Review docs
• Post review

Note: I could not execute Python in this environment (perceval isn't installed and python invocations are blocked), so everything below is static analysis — no test run. CI "Docs Validation" is green; "CI Tests" was still in progress.

---

The core idea is sound and the implementation is compact and well documented. The global-phase folding is mathematically correct: subtracting alpha from each of the m phases of a full phase layer multiplies the block unitary by e^{-i·alpha}, which is exactly the required correction, and the trace-based alpha is the least-squares optimal estimate. Keeping the fast path when no phase noise is configured, and excluding PERM, are the right calls.

The main problems are in the new tests and in the accuracy/determinism of the numerical fit.

High

1. test_reservoir_classifier_with_phase_noise_on_gpu asserts behavior that the cache makes impossible — tests/models/test_reservoir_classifier.py

logits1 = model.predict(X)
logits2 = model.predict(X)
assert not torch.allclose(logits1, logits2, rtol=0.1)

ReservoirClassifier defaults to cache=True (merlin/models/reservoir_classifier.py:120), and predict(X) on the same X that was passed to fit_reservoir takes the cache-hit branch at reservoir_classifier.py:731-735, returning self._fit_quantum_cache.clone(). Both calls return the identical tensor, so not torch.allclose(...) is False and the assert fails. It also directly contradicts test_reservoir_classifier_with_phase_imprecision_and_error_combined in the same file, which asserts torch.allclose(logits1, logits2) for the same noise configuration. The "Tests pass on GPU" box is unchecked, so this was never executed. Separately, rtol=0.1 inside a not allclose is a very weak difference check even if the cache were disabled.

To actually exercise per-call stochastic sampling, build the model with cache=False (or call model._quantum_layer(...) directly on fresh inputs).

Fix this →

2. The three new CPU reservoir tests pass for the wrong reason, and their docstrings state incorrect behavior — tests/models/test_reservoir_classifier.py

test_reservoir_classifier_with_phase_error_only documents:

Without apply_phase_error=True, stochastic sampling is NOT triggered … For typical use (fit/predict without explicit stochastic sampling), only phase_imprecision (if configured) affects outputs.

That is not what the code does. ComputationProcess.compute() routes to _compute_phase_error_probabilities() whenever has_phase_error(noise_groups) is true (merlin/core/process.py:978-979), and that method calls to_tensor(..., apply_phase_error=True) in a loop over n_phase_error_samples (process.py:881-894). So a NoiseModel(phase_error=0.1) does sample stochastically on every uncached forward. The tests observe identical results only because of the feature cache (see #1). As written, these three tests assert cache behavior, not noise behavior, and they encode a wrong mental model in the docstrings — which is worse than no test, since a future regression in the noise path would still pass. Per AGENTS.md ("keep tests close to the behavior they protect"), these should either use cache=False and assert variability, or be honest that they are cache tests.

3. The decomposition introduces an uncontrolled ~1e-3 noise floor that can exceed the noise being modelled — merlin/pcvl_pytorch/locirc_to_tensor.py:149

CircuitOptimizer() is used with default settings, giving fidelity error ~1e-6, i.e. matrix entries off by ~1e-3 (as your own docs state). For a user modelling realistic hardware imprecision — e.g. phase_imprecision = 2π/1024 ≈ 6e-3 rad, or phase_error of a few mrad — the fit error is the same order as, or larger than, the physical effect being simulated. The "noisy" result is then dominated by an artifact of the fit. Two options worth considering:

• pass a tighter threshold to CircuitOptimizer(...), at the cost of build time; or
• use the exact algebraic Clements decomposition (pcvl.Circuit.decomposition(..., shape=InterferometerShape.RECTANGLE)), which reproduces the target to numerical precision and does not depend on optimizer convergence.

At minimum, warn when the achieved fit error is comparable to the configured phase_imprecision/phase_error.

4. No post-fit verification; alpha is only meaningful if the fit succeeded — locirc_to_tensor.py:148-161

The only failure detection is except RuntimeError around optimize_rectangle. Please confirm that this perceval version raises rather than logging a warning and returning a poor fit — if it warns, a bad decomposition silently replaces the user's unitary, and alpha = angle(trace(target^H @ fitted)) becomes meaningless (the trace can be near zero, making the phase arbitrary). AGENTS.md is explicit here ("Do not hide failures"). A three-line guard makes this robust regardless of perceval's behavior:

fitted = np.asarray(mesh.compute_unitary(), dtype=complex)
overlap = np.trace(target.conj().T @ fitted) / component.m
if abs(overlap) < 1 - _FIT_TOLERANCE:
    raise RuntimeError(...)
alpha = float(np.angle(overlap))

Medium

5. Non-determinism across constructions is understated. optimize_rectangle starts from random points drawn from the global NumPy RNG, so (a) each CircuitConverter build produces different mesh phases, (b) torch.manual_seed does not control it, and (c) building a layer perturbs the user's global NumPy RNG state. The docs say "the mesh phases differ between constructions … but the reproduced unitary does not" — that holds only up to the ~1e-3 fit error, and it is false once phase_imprecision quantization is applied, since quantization is a discontinuous function of the individual mesh phases. Practical consequence: save/reload or rebuild of the same model yields a different noisy unitary. Consider a module-level cache keyed on the matrix bytes, and/or an explicit seed, and correct the doc wording.

6. No memoization of the decomposition. Every CircuitConverter construction re-runs the optimizer for every Unitary, including identical matrices, and ReservoirClassifier._rebuild_quantum_layer() re-runs it whenever layer.noise is assigned. With the quoted ~0.2 s for 10 modes and superlinear growth, a 16–20 mode reservoir rebuilt a few times becomes noticeable. A dict keyed on (matrix.tobytes(), m) would be cheap and would also mitigate #5.

7. Only PS phases receive noise — the mesh may hide phases inside BS. _compute_tensor(BS) (locirc_to_tensor.py:782-853) applies neither quantization nor perturbation; any phase carried by a BS parameter (theta, phi_tl, …) in the template used by optimize_rectangle stays perfectly noiseless. The tests only assert any(isinstance(c, PS) for _, c in converter.list_rct), which cannot distinguish "all mesh phases are noisy" from "only the output layer is noisy". Since the stated goal is "phase noise on every physical phase", please add a test asserting the expected count of noise-carrying PS for an m-mode mesh (e.g. 2·m(m−1)/2 + m for the phase-last MZI template). That test also guards the p.name.startswith("phL") prefix coupling at line 163, which is tied to perceval's template naming (the dependency is only pinned as perceval-quandela>=1.2.1).

8. Unhandled edge cases. A 1-mode non-polarized Unitary (a bare phase) goes straight into optimize_rectangle, which has no MZI to place — likely an obscure crash rather than the clear error the guard at line 164 intends. A trivial component.m == 1 branch emitting a PS would cover it. Large-m build cost also deserves either a documented scaling note or a warning above some mode count.

Low / nits

• Trailing whitespace: tests/models/test_reservoir_classifier.py has 40 lines with trailing spaces, which ruff format would strip and W293 would flag. CI only lints merlin/ (.github/workflows/ci.yml:56-57), so this won't fail the build, but the "Code formatted" checkbox isn't accurate for this file.
• Arbitrary assertions: assert feature_diff < 1.0 and assert feature_diff > 0.0 don't encode a meaningful property. Also, the new reservoir tests don't seed NumPy, so both the random reservoir unitary and the fit vary run to run — feature_diff < 1.0 is a flakiness candidate.
• Tolerance mismatch: the new test module documents "entries agreeing to about 1e-3" but sets _ELEMENTWISE_ATOL = 1e-2.
• Docs timing claim: "about 0.2 s for a 10-mode unitary" is machine-specific; a scaling statement would age better.
• _decompose_unitaries drops circuit.name on the rebuilt circuit (cosmetic; affects display only).
• Behavior change worth a release note: existing users with only phase_imprecision set and a Unitary in the circuit will now get different numerical results and a slower build.

What's good

• Fast path preserved and explicitly regression-tested (test_reservoir_unitary_layer_neutral_noise_matches_no_noise, test_no_noise_converter_keeps_unitary_fast_path).
• PERM exclusion, with a test asserting object identity.
• Global-phase compensation for sub-mode blocks is correct and has a dedicated test.
• Error messages name the component and its modes; polarized and non-convergence paths both tested.
• Gradient flow through input parameters checked end-to-end on the decomposed circuit.
• Docs in noisy_simulations.rst explain the "why" (hardware phases inside black boxes) clearly.
  • Branch PML-361-circuit-noise-from-unitary

[CassNot]

@claude review

[claude]

Claude finished @CassNot's task in 5m 47s —— View job

---

Code review — PML-361 (phase errors on Unitary-based circuits)

• Gather context (AGENTS.md, PR diff vs origin/release/0.4.1)
• Review _decompose_unitaries() + integration in CircuitConverter
• Re-check the fixes from 33a280a, 0f7c64f, cfee3fa against my previous review
• Review tests and docs
• Post review

Static analysis only — perceval isn't installed here and python/ruff invocations are blocked in this sandbox, so nothing below was executed. CI: "Docs Validation" green, "CI Tests" still in progress at the time of writing.

---

Good news first: the three findings I flagged last round are genuinely fixed. The reservoir tests now use cache=False, so they exercise the real per-call sampling path (_compute_quantum_features line 740-744 falls through to _encode_quantum on every call), and the misleading docstrings are gone. The post-fit overlap guard (locirc_to_tensor.py:271-279) and the memoization (_DECOMPOSITION_CACHE) both landed, and the m == 1 and m >= 16 edge cases are handled.

The dominant problem now is that the last commit — cfee3fa "no phase error imprecision on BS" — is an incomplete revert. It removed the BS call site but left the machinery and all the documentation behind, so the code and its own docs now describe opposite behavior.

High

1. _apply_phase_noise is dead code — merlin/pcvl_pytorch/locirc_to_tensor.py:922-968

cfee3fa removed its only call site (the BS loop). grep -rn "_apply_phase_noise" merlin/ returns exactly one hit: the def. The PS overload never migrated to it — it still carries the inline copy of the same quantization + perturbation logic at lines 1148-1176. So there are now 46 lines of unused code duplicating live logic, and its docstring claims something false:

This helper is used by the :meth:`_compute_tensor` overload for ``PS``.

Either delete the helper, or (better, per AGENTS.md "Do not duplicate logic") have the PS overload actually call it.

Fix this →

2. The BS clause in is_phase_error_sensitive is now a pure performance regression — locirc_to_tensor.py:718-722

is_phase_error_sensitive = (
    isinstance(c, PS) and (...)
) or (isinstance(c, BS) and self._phase_error > 0.0)

With that clause, every constant BS is kept as a live component instead of a precomputed tensor whenever phase_error > 0 — but _compute_tensor(BS) (lines 1010-1070) now applies no noise at all, so the extra work changes nothing numerically. Two costs, both on the exact path this PR adds:

• every BS tensor is rebuilt on every to_tensor() call, and _compute_phase_error_probabilities calls to_tensor(apply_phase_error=True) n_phase_error_samples times per forward (merlin/core/process.py:881);
• worse, the second pass in _convert_circuit (lines 730-760) can only fuse adjacent numeric tensors. Keeping the BS as components blocks their modes and defeats the merge for the whole mesh. A decomposed 10-mode unitary has 45 MZIs ≈ 90 BS, all of them now unmergeable, where before the whole Unitary was a single constant tensor.

The accompanying comment also states the opposite of what the code does:

# BS: all five BS parameters (theta, phi_tl, phi_bl, phi_tr, phi_br)
#     are physical phases and receive the same global phase_error, so
#     a constant BS must also stay dynamic when phase_error is active.

Drop the isinstance(c, BS) clause and the comment.

Fix this →

3. Four places still document BS noise that no longer exists. After #2 the behavior is "PS only", but these all still claim otherwise:

Location	Stale claim
locirc_to_tensor.py:141-143 (_decompose_unitaries docstring)	"reaches every physical phase … including the theta and phi parameters of every BS in the mesh"
locirc_to_tensor.py:718-721 (comment)	"a constant BS must also stay dynamic"
docs/source/quantum_expert_area/noisy_simulations.rst:112	"so the configured phase noise applies to every physical phase"
tests/pcvl_pytorch/test_unitary_decomposition_noise.py:253	"phase_imprecision=0.1 quantizes every BS and PS phase in the mesh"

The rst line matters most — it's the user-facing promise, and it's the PR's stated goal. If the MZI template used by optimize_rectangle keeps its two tunable phases as PS (which the phL* output-layer naming suggests), then "PS only" is complete in practice and the docs just need rewording to say so explicitly. If any tunable phase lives in a BS parameter, then it now receives no noise and the goal is only partly met. Which is the case is worth stating in the docs either way, since it's the difference between "all mesh phases are noisy" and "only the output layer is noisy".

4. Neither new warning has a test, and neither new failure branch does. AGENTS.md is explicit: "If new warning behavior is added, add or update tests that assert the warning is emitted." grep over both test files finds no pytest.warns for:

• the m >= 16 slow-decomposition warning (line 187);
• the fit-error-vs-noise-scale warning (line 289).

Also untested: the overlap-quality RuntimeError (line 272), the phL count-mismatch RuntimeError (line 314), and the new component.m == 1 branch (line 170) — test_polarized_unitary_raises_value_error builds a 1-mode circuit but raises on polarization before reaching it, so the bare-phase path has zero coverage.

Medium

5. The fit-error warning is mis-scaled by ~3 orders of magnitude and will essentially never fire — locirc_to_tensor.py:281-299

fit_error = 1.0 - abs(overlap)
if fit_error > 0.1 * noise_scale:

1 - |overlap| is second order in the phase discrepancy, while noise_scale is first order. For per-entry error σ, |overlap| ≈ 1 - σ²/2. Your own docs state entries agree to ~1e-3, i.e. fit_error ≈ 5e-7. The condition then needs noise_scale < 5e-6 rad to trigger — but the artifact it's meant to catch (σ ≈ 1e-3) already dominates a phase_error of 1e-4 rad, silently. Compare like with like:

entry_error = float(np.sqrt(2.0 * max(0.0, 1.0 - abs(overlap))))  # ≈ per-phase discrepancy, rad
if entry_error > 0.1 * noise_scale:

Same dimensional issue in _FIT_TOLERANCE = 1e-3: as a bound on 1 - |overlap| it admits σ ≈ 0.045 rad, which is a much looser gate than the docstring's "overlap should be > 0.999" implies (with fidelity error ~1e-6 you should see > 0.999999). The gate is still a real improvement over no check — it's the warning that's ineffective.

Fix this →

6. _DECOMPOSITION_CACHE is an unbounded process-global with no way to clear it — locirc_to_tensor.py:47

The memoization is the right fix for the RNG/determinism problem, but as written:

• No eviction. Every distinct matrix permanently retains a full mesh Circuit (≈ m(m−1) components; ~380 for m=20). A ReservoirClassifier sweep or ensemble — random reservoir unitary per trial, which is the headline use case — leaks one entry per trial for the process lifetime. A bounded LRU (collections.OrderedDict + maxsize, or functools.lru_cache over a small keyed helper) plus a documented clear_decomposition_cache() would fix both this and test isolation.
• The determinism claim overreaches. The docstring says "(b) rebuilding or reloading a model always produces the same mesh phases". That holds within a process; a save + reload in a fresh interpreter re-runs optimize_rectangle from a fresh global NumPy RNG and gets different phases — and since phase_imprecision quantization is discontinuous in those phases, a different noisy unitary. Please scope the claim to a single process. Truly reproducible results need either an explicit seed or the exact algebraic Clements decomposition (pcvl.Circuit.decomposition(..., shape=InterferometerShape.RECTANGLE)), which would also remove the ~1e-3 artifact behind move tests to root, add license header in all files #5.
• No test. The cache is the mechanism behind the determinism guarantee and has no direct coverage. It also has no isolation fixture, so test_non_convergence_error_names_component (which monkeypatches optimize_rectangle to raise, and doesn't seed NumPy) depends on its unseeded matrix never being a cache hit. An autouse fixture clearing the dict would make that robust. Note that test_circuit_converter_with_phase_noise_on_gpu's final atol=1e-12 CPU-vs-GPU comparison silently relies on the cache — worth a comment, since it looks like a device-parity assertion but is really a cache-determinism assertion.

7. The m >= 16 warning fires before the cache lookup — locirc_to_tensor.py:185-193

It sits above the cache_key in _DECOMPOSITION_CACHE check, so a cached 16-mode unitary still warns "this may take several seconds … The result is cached after the first build" on every subsequent build, when the build is in fact instant. Move it below the cache check. (The fit-error warning has the mirror-image asymmetry — it only fires on a cache miss, so a genuine accuracy problem is reported once and then goes quiet for the rest of the process. Worth a deliberate decision on which of the two you want.)

Low / nits

• ruff format was still not run on tests/models/test_reservoir_classifier.py — 41 lines with trailing whitespace (grep -c '[ \t]$'). CI only lints merlin/, so it won't fail the build, but the "Code formatted" checkbox isn't accurate for that file.
• Cosmetic churn in _compute_tensor(BS). After cfee3fa the raw = ...; param_values.append(raw) indirection (lines 1013-1026) is a no-op leftover, and the reformatted unitary_tensor = (unitary_tensor .unsqueeze(0) ...) chain at lines 1060-1065 is a leading-dot style ruff format would collapse. Reverting both keeps the diff to the actual change (AGENTS.md: "Do not rewrite unrelated code").
• assert feature_diff < 1.0 / assert feature_diff > 0.0 (test_reservoir_classifier.py) still don't encode a property; the second is trivially true for any float comparison.
• not torch.allclose(a, b, rtol=0.1) — this passes largely because standardized features contain near-zero entries, where 0.1 * |b| collapses to atol. It reads as a strong check but is really an incidental one. (a - b).abs().max() > 1e-6 would state the intent (fresh noise sample per call) directly and won't drift.
• Tolerance mismatch persists: test_unitary_decomposition_noise.py documents "entries agreeing to ~1e-3" but sets _ELEMENTWISE_ATOL = 1e-2.
• phL prefix coupling (line 306) is still tied to perceval's internal template naming with only perceval-quandela>=1.2.1 pinned. The count-mismatch RuntimeError guards it, which is the right shape — just make sure the pin reflects the version that guarantees the naming.
• _decompose_unitaries drops circuit.name on the rebuilt circuit (display only).
• Release note still worth adding: existing users with only phase_imprecision set and a Unitary in their circuit get different numerical results and a slower first build.

What's good

• The cache=False rework of the reservoir tests is exactly right — they now exercise the real stochastic path instead of the feature cache.
• The overlap-quality gate makes a bad fit loud instead of producing a meaningless alpha, per the failure policy.
• Global-phase folding for sub-mode blocks is correct and has a dedicated test; PERM exclusion is tested by object identity.
• Fast path preserved and regression-tested from both ends (test_reservoir_unitary_layer_neutral_noise_matches_no_noise, test_no_noise_converter_keeps_unitary_fast_path).
• Memoization is the right call and the cache docstring reasons carefully about thread safety and RNG perturbation.
• Gradient flow through input parameters verified end-to-end on the decomposed circuit, plus unitarity checks.
  • Branch PML-361-circuit-noise-from-unitary

[CassNot]

I agreed with all Claude's comments (High) and addressed them. I also added more thorough GPU tests to ensure compatibility and they all run.
@claude review

[claude]

Claude finished @CassNot's task in 6m 50s —— View job

---

Code review — PML-361 (phase errors on Unitary-based circuits)

• Gather context (AGENTS.md, pytest.ini, PR diff vs origin/release/0.4.1)
• Re-check previous High findings against f7cd02b, 1e7a71e, 6e71c7c, 096f530, 7ad6beb
• Review _decompose_unitaries() + integration
• Review tests (incl. new GPU tests) and docs
• Post review

Static analysis only — perceval isn't installed here and python / pip invocations are blocked in this sandbox, so nothing below was executed. Where I claim a perceval API behavior I cite in-repo evidence rather than memory.

---

The three previous High findings are genuinely fixed: _apply_phase_noise is now actually called from the PS overload (locirc_to_tensor.py:1197), the isinstance(c, BS) clause is gone from is_phase_error_sensitive (line 761-763) so the numeric-merge pass works again, and every new warning and failure branch now has a test. The phase-aligned Frobenius residual (line 331-335) is exactly the right first-order measure — my earlier scaling complaint is fully addressed.

test_decomposed_mesh_keeps_all_phase_shifters_noise_sensitive also settles the question I left open last round: 2·m(m−1)/2 + m PS for m=4 means every tunable mesh phase is a PS, so "PS only" is complete coverage. Which makes it all the more unfortunate that the docs now say the opposite (see #2).

Two High items remain.

High

1. The memoization never hits — every build still re-runs the optimizer, and the entry is deleted each time — merlin/pcvl_pytorch/locirc_to_tensor.py:272-284

Line 357-361 fixes every mesh parameter via parameter.fix_value(value). Then on the next build, the cache-hit branch does:

cached_parameters = cached_mesh.get_parameters(all_params=False)   # line 272
cached_output_phases = [p for p in cached_parameters if p.name.startswith("phL")]
if len(cached_output_phases) == component.m:
    new_circuit.add(r[0], cached_mesh, merge=True); continue
del _DECOMPOSITION_CACHE[cache_key]                                 # line 284

get_parameters(all_params=False) returns only non-fixed parameters, so on a fully-fix_valued mesh it returns []. Two independent proofs in this repo:

• benchmarks/benchmark_unitary_conversion.py:94-96 — for param in circ.get_parameters(): param.fix_value(...) followed by assert len(circ.get_parameters()) == 0.
• This PR's own tests/pcvl_pytorch/test_unitary_decomposition_noise.py:85 — assert decomposed.get_parameters() == [] on the decomposed circuit.

So cached_output_phases is always [], len([]) == component.m is always false for m ≥ 2, the entry is deleted, and the optimizer re-runs. The cache is a permanent miss that also churns its own storage. Consequences, in order of how much they matter:

• All three guarantees in the docstring (lines 158-162) are false. The global NumPy RNG is perturbed on every build, not once per matrix; rebuilds do not produce the same mesh phases; phase_imprecision quantization is not consistent across builds. This was the whole point of adding the cache last round.
• ReservoirClassifier loses its intended left/right symmetry under noise. reservoir_classifier.py:367-368 builds pcvl.Unitary(matrix) and pcvl.Unitary(matrix.copy()) — byte-identical matrices, same cache key. With a working cache both get the same mesh; as-is they get two independent fits, so under phase_imprecision the "symmetric output interferometer" quantizes to a different unitary than the input one. _rebuild_quantum_layer() re-runs both optimizations on every layer.noise assignment.
• test_circuit_converter_with_phase_noise_on_gpu:461-463 is now a coin flip, not a device check. It compares two independently-fitted meshes, each ~1e-3 per-entry off the target, at atol=1e-3. That is why 7ad6beb had to loosen it from 1e-12 — the loosening was treating the symptom. With the cache working the two would be bit-identical and 1e-12 would hold.

Fix is one word (all_params=True at line 272), plus a test — the cache has zero direct coverage today, which is why this shipped. A test asserting that a second _decompose_unitaries call on the same matrix yields bit-identical phL values (and, with a counting monkeypatch, does not call optimize_rectangle twice) would have caught it.

Fix this →

2. Three places now document the exact opposite of what the code does — mesh PS do receive noise

Location	Claim	Reality
locirc_to_tensor.py:142-143	"the fixed mesh parameters themselves are not tunable and receive no phase noise"	wrong
docs/source/quantum_expert_area/noisy_simulations.rst:113-115	"This makes the unitary's structure explicit for analysis. Phase noise applies only to explicitly defined tunable pcvl.PS components in the circuit"	wrong
tests/pcvl_pytorch/test_unitary_decomposition_noise.py:413	"Decomposed mesh parameters are fixed constants and receive no phase noise"	wrong

_compute_tensor(PS) reaches _apply_phase_noise on both branches — the is_variable one at line 1168 and the fixed-value _value one at line 1181-1184 — quantization and perturbation are applied unconditionally at line 1197. The PR's own tests confirm it: line 141-160 counts 16 noise-carrying PS in a 4-mode mesh, and line 178-188 asserts the quantized mesh output deviates from the target. The class docstring at locirc_to_tensor.py:459-466 says the correct thing ("Phase noise is then applied to the PS components of the mesh") and directly contradicts the function docstring 320 lines above it.

The rst line is the one that actually matters: it is the user-facing promise, and as written it tells a reader that this PR's feature does not do what the PR title says. Since #1 confirms all mesh phases are PS, the honest wording is "phase noise reaches every tunable phase of the mesh, which are all pcvl.PS; pcvl.BS parameters carry no noise, and the mesh's beam splitters are fixed 50:50 couplers with no tunable phase".

Fix this →

Medium

3. The cache write is not guarded by optimizer_is_default, so a patched optimizer poisons a process-global. The read is guarded (line 251), the write at line 366 is not. Concretely: test_fit_error_vs_noise_scale_warning warns rather than raises, so its MockMesh is deep-copied into _DECOMPOSITION_CACHE under the seed-4 matrix key — the same key used by test_phase_error_sampling_varies_and_is_reproducible, which also does np.random.seed(4) + _single_unitary_circuit(4). File order saves it today (line 163 runs before line 275); any reordering (-k selection, split runs, random ordering) serves a mock mesh to a real test, and the mock's get_parameters returns 4 phL entries so the reuse validation would pass it. Seeds 5 and 6 collide the same way but those tests raise before the write. Guard the write with optimizer_is_default, and add an autouse fixture clearing the cache — the same fixture also removes the current unseeded dependence in test_non_convergence_error_names_component.

4. SUPPORTED_COMPONENTS lost its docstring — locirc_to_tensor.py:46-83. The new _DECOMPOSITION_CACHE (48-71) and _DEFAULT_OPTIMIZE_RECTANGLE (73-74) assignments-plus-docstrings were inserted between SUPPORTED_COMPONENTS = (...) and its docstring, so the "Tuple of quantum components supported by CircuitConverter" block at lines 75-83 is now a stray literal trailing _DEFAULT_OPTIMIZE_RECTANGLE's own docstring — documenting nothing, and invisible to autodoc for either symbol. Move the new globals below line 83.

5. The m >= 16 warning still fires above the cache lookup — lines 230-238 precede the cache_key in _DECOMPOSITION_CACHE check at 251. Masked today by #1, but once the cache works a cached 16-mode unitary will warn "may take several seconds … The result is cached after the first build" on every rebuild, when the build is instant. The residual warning has the mirror asymmetry (miss-only, so a genuine accuracy problem is reported once per process then goes quiet) — worth a deliberate decision on which side each belongs. Separately, test_slow_decomposition_warning_for_large_unitary runs a real 16-mode optimization, which by the warning's own text is "several seconds" of suite time; monkeypatching the optimizer keeps the assertion and drops the cost.

6. _FIT_TOLERANCE is still a second-order gate behind a now-first-order warning — lines 85-94, 304. The hard RuntimeError compares 1 - |overlap| against 1e-3; since |overlap| ≈ 1 - σ²/2, that admits a per-entry discrepancy of σ ≈ 0.045 rad. The warning at line 336 now correctly uses the first-order phase-aligned residual. Net effect: a band of genuinely bad fits warns but never errors. Also _FIT_TOLERANCE's docstring says "the overlap should be > 0.999" where the real figure is > 0.999999 — the stated bound is 3 orders looser than the behavior it describes.

7. _DECOMPOSITION_CACHE is still unbounded with no clear API (carried over). Once #1 is fixed this becomes load-bearing: one full mesh Circuit (≈ m(m−1) components) retained per unique matrix for the process lifetime. A ReservoirClassifier sweep or ensemble draws a fresh random unitary per trial — one leaked entry per trial. A bounded OrderedDict LRU plus a documented clear_decomposition_cache() covers both this and test isolation (#3). Also please scope the "rebuilding or reloading a model always produces the same mesh phases" claim (line 159-160) to a single process — a save/reload in a fresh interpreter re-fits from a fresh RNG.

Low / nits

• ruff format still not run on the touched test files. tests/models/test_reservoir_classifier.py has 41 trailing-whitespace lines, and the new tests/pcvl_pytorch/test_unitary_decomposition_noise.py isn't formatted either (lines 153-155 collapse to one line, the 267-268 circuit.add(...) split, over-length monkeypatch.setattr at 311/343/373, and 375). CI only lints merlin/ (.github/workflows/ci.yml:56-57), so the build won't fail — but the "Code formatted (ruff format)" checkbox isn't accurate.
• test_one_mode_unitary_creates_single_ps:249 assigns phase_value = np.pi / 4 and never uses it (ruff F841). More to the point, the test asserts only that a PS is emitted, never that its phase equals angle(target[0, 0]) — which is the part of line 224 worth protecting.
• _apply_phase_noise docstring first line: "Apply configured phase noise to a phase tensor in-place (returns new tensor)" — self-contradictory; nothing is in-place.
• assert feature_diff > 0.0 (combined reservoir test) is still trivially true for two independently-sampled noisy runs.
• not torch.allclose(a, b, rtol=0.1) is used in five places as the "differs" check. On standardized features 0.1 · |b| collapses to atol near zero entries, so it reads much stronger than it is. (a - b).abs().max() > 1e-6 states "fresh noise sample per call" directly.
• Tolerance mismatch persists: the new module docstring says "entries agreeing to ~1e-3" while _ELEMENTWISE_ATOL = 1e-2.
• pytest.ini:5-6 sets filterwarnings = error, so both new UserWarnings are exceptions under test. The reservoir tests run at m = 3 (reservoir_classifier.py:197, encoded_input_features + 1 with PCA(n_components=2)), residual ≈ 1.7e-3 against thresholds of 5e-3–1e-2, so there's headroom — but a future test at realistic mrad-scale phase_error would hard-fail rather than warn. Worth knowing this is the chosen behavior.
• The alpha fold at line 357-361 can push a phL value negative and relies on perceval's Parameter being periodic so fix_value wraps rather than range-failing. Worth one inline sentence, since it's an implicit dependency on the template's parameter construction — same category as the phL prefix coupling (guarded by the count check, good) with only perceval-quandela>=1.2.1 pinned.
• _decompose_unitaries still drops circuit.name on the rebuilt circuit (display only).
• Release note still worth adding: existing users with only phase_imprecision set and a Unitary in their circuit get different numerical results and a slower build.

What's good

• _apply_phase_noise is now genuinely shared by the PS overload — the dead code and the duplicated quantization/perturbation logic are both gone.
• Dropping the BS clause from is_phase_error_sensitive restores the adjacent-numeric-tensor merge for the whole decomposed mesh; that was the biggest hidden cost.
• The phase-aligned Frobenius residual is dimensionally correct and detects entry-wise error the trace overlap can cancel.
• Every new warning and failure branch is tested: m >= 16, fit residual, overlap gate, phL count mismatch, m == 1, polarized, non-convergence. That's the full set AGENTS.md asks for.
• test_decomposed_mesh_keeps_all_phase_shifters_noise_sensitive with the exact 2·m(m−1)/2 + m count is precisely the right test, and it also answers the coverage question definitively.
• Read-side cache bypass under a replaced optimizer, and cached-mesh unitarity revalidation, are both careful defensive additions — the reuse validation is the right shape, it just checks the wrong parameter set (Update issue templates #1).
• Global-phase folding for sub-mode blocks remains correct and dedicated-tested; PERM exclusion tested by object identity; fast path regression-guarded from both ends.
  • Branch PML-361-circuit-noise-from-unitary

[LF-Vigneux]

Everything seems good