cirq.unitary(…) can return a non-unitary in some cases

[mhucka]

@pavoljuhas noted in a comment on PR quantumlib/qsim#857 that as things stand now, cirq.unitary(…) may return a non-unitary value, which means that it's possible for cirq_unitary(cirq.unitary(val)) to fail when cirq.unitary(val) works.

This became possible after PR #7419 added checks in the case of NumPy arrays. This is arguably the right behavior, in the sense that cirq.unitary(…) is supposed to return a unitary but previously would simply return the argument unverified when the argument was a NumPy array. Now it won't do that, but the other strategies used by cirq.unitary(…) don't do similar checks, and that leads to the possibility of new inconsistent behaviors.

[daxfohl]

PR #7419 added checks in the case of NumPy arrays

Kind of: that PR added the ability to call cirq.unitary on numpy arrays at all. Prior to that PR, calling cirq.unitary on a raw numpy array would throw regardless of whether it was unitary or not. Now it succeeds for raw unitary arrays and only throws for non-unitaries. It's a ValueError now for non-unitary arrays, versus pre-7419 it would throw a TypeError, but otherwise should be the same. I don't see anywhere in the test that checks the type of error, so there's something else going on that is causing the test to crash.

[daxfohl]

Oh I see: if you go back and look at the failing call stack https://github.com/quantumlib/qsim/actions/runs/16846849990/job/47727587531#step:6:1003, you see it hits this line 97 from mixture_protocol.py, that was also changed in 7419: https://github.com/quantumlib/Cirq/pull/7419/files#diff-39c5258dee971ba4f4af1c6734a57cc85d175c824b4b11f128623065677f8200.

You can see that previously, it'd return a raw array as-is but call unitary for anything else. Since the PR allowed unitary to handle raw arrays, that distinction was deemed no longer necessary, so the special case was removed. But of course that assumed the mixture components' arrays are actually unitary. Mixture components should all be unitary, so I think the test itself is probably malformed? But if there's a need to support mixtures with non-unitary components, then reverting line 97 there would fix the problem. (Or of course, removing the unitarity check from unitary would fix it too).

I'd say my first choice would be to leave this as-is: the test case in qsim uses a malformed mixture, and so the fix that was applied to the qsim test seems sufficient. Second choice would be to have cirq.unitary accept an optional tolerance parameter, and only check unitarity if that tolerance parameter is provided (the downside being that then you feel like you need to provide a tolerance parameter everywhere, for consistency, which, can of worms). Third choice would be to remove the unitarity check from cirq.unitary when the input is a raw array; it's probably fine, but seems like it would be surprising from an end-user perspective. Fourth choice would be to undo #7419 altogether; it'd fix the problem, but the functionality in that PR is useful. Last choice would be to have it check unitarity on all return values; it's expected that the implementations of _unitary_ etc are indeed unitary, and checking those for unitarity before returning it would be prohibitively expensive.

[mhucka]

Thanks for writing a succinct summary of the problem! Yes, the problem was indeed due to combination of changes (to unitary and to mixture protocol).

I agree with your first choice. It really seemed like the qsim test took a shortcut in the implementation of the simple mixture object – not surprising, since it's only a test case.

[pavoljuhas]

I think it is better to go for option 3, ie, removing the is_unitary check here

Cirq/cirq-core/cirq/protocols/unitary_protocol.py

Lines 117 to 118 in dc69fa1

	if isinstance(val, np.ndarray):
	if not linalg.is_unitary(val):

This would resolve the inconsistency of unitary(val) working and unitary(unitary(val)) failing. There are several places in the cirq test suite that have cirq.unitary return a non-unitary matrices so adding extra checks per options 2 or 5 would break them and would introduce new incompatible behavior. The docstring of cirq.unitary states the value of obj._unitary_() is returned unchecked -

Cirq/cirq-core/cirq/protocols/unitary_protocol.py

Lines 89 to 90 in dc69fa1

	- The value has a `_unitary_` method that returns something besides None or
	NotImplemented. The matrix is whatever the method returned.

In addition cirq.unitary has the default argument which is returned unchecked and can be anything. My reading of cirq.unitary docs is it does the best effort to compute a unitary, but does not guarantee that value to be such; the is_unitary is a proper check for that.

PS: Here is a patch to test-add is_unitary check

diff --git a/cirq-core/cirq/protocols/unitary_protocol.py b/cirq-core/cirq/protocols/unitary_protocol.py
index 0e779472..1f5254e3 100644
--- a/cirq-core/cirq/protocols/unitary_protocol.py
+++ b/cirq-core/cirq/protocols/unitary_protocol.py
@@ -131,2 +131,4 @@ def unitary(
         if result is not NotImplemented:
+            if not linalg.is_unitary(result):
+                raise ValueError("The provided result is not unitary.")
             return result

and here are tests that it would break:

cirq-core/cirq/testing/circuit_compare_test.py::test_assert_has_consistent_apply_unitary
cirq-core/cirq/protocols/pauli_expansion_protocol_test.py::test_pauli_expansion[val3-expected_expansion3]
cirq-core/cirq/protocols/pauli_expansion_protocol_test.py::test_pauli_expansion[val2-expected_expansion2]

[daxfohl]

I think it's fine, though one extra thing to consider: #7537 added a tolerance check to the kraus/channel protocol. I didn't notice that until yesterday. But in that light, do you think it would make option 2 more favorable for consistency sake? I wonder if there are similar rounding issues in the "unitary_from_apply_unitary", that just haven't been noticed yet.

But otherwise I think option 3 is fine too.

[daxfohl]

Nvm the above. Looking at that PR, the atol is used to filter the "rounding error from zero" terms from the calculated Kraus representation, so serves a different purpose than option 2 here. If both protocols have an atol but they're used for different things, that makes it even more confusing to use. So option 3 seems good.

Plus I like that it allows unitary(unitary(x)) == unitary(x), which seems like a reasonable thing to expect from most protocols.