Fix bogoliubov_transform for spin-preserving Gaussian transformations

[blazingphoenix7]

Fixes #776.

bogoliubov_transform checks whether the transformation matrix is spin-block-diagonal before it ever looks at the matrix shape, and that check only makes sense for a square N x N matrix. For the N x 2N matrix you get out of a general (non-particle-conserving) Bogoliubov transformation it ends up inspecting the wrong blocks, decides that a transform which doesn't actually mix spin is "block diagonal", and then slices it into pieces of the wrong shape. The matrix from the issue dies with:

ValueError: Bad shape for transformation_matrix. Expected (1, 1) or (1, 2) but got (1, 3).

What I changed

I moved the block detection into a helper, _spin_blocks, that handles both shapes. For the N x 2N case, write the matrix as W = (W1 W2) (it acts on all the creation operators followed by all the annihilation operators). The transform leaves the two spin sectors alone exactly when both W1 and W2 are block diagonal, and in that case each sector is its own (N/2) x N Gaussian transform built from the matching diagonal blocks of W1 and W2. The recursion that was already there handles the rest. As a bonus, for spin-symmetric systems (BCS-style pairing Hamiltonians, say) you now do two half-size Givens decompositions instead of one full-size one.

The part that took some digging

I didn't want to call this correct based on eigenstate energies alone, so I checked the circuit at the operator level: does it actually give U a_p† U⁻¹ = b_p† for every mode? That caught something the energy checks can't.

Under Jordan-Wigner a spin-down ladder operator drags a Z string across all of the spin-up qubits, so it anticommutes with the spin-up parity operator. If the spin-up circuit flips spin-up parity (which happens whenever its Gaussian decomposition uses an odd number of particle-hole transforms), then splitting the sectors apart flips the sign of every spin-down operator. What makes it sneaky is that it's invisible to energy tests: starting from a fixed occupation state, that bad sign is just a global phase, so the state still comes out as the correct eigenstate with the correct energy even though the underlying unitary is wrong.

So when the spin-up circuit flips parity I add a layer of Z gates on the spin-down qubits to put the signs back (_preserves_parity decides whether that's needed). With a fixed computational-basis initial state it's only a global phase, so I skip it there.

I also ran into the general Gaussian decomposition not reliably diagonalizing exactly-block-structured matrices (U† H U comes out non-diagonal for some perfectly valid block inputs), so I went with always splitting these instead of trying to detect when it's safe to fall back. Splitting is the more correct path here, not just the faster one.

Tests

• test_bogoliubov_transform_spin_block_gaussian_regression: the actual matrix from The implementation of bogoliubov for a general fermionic gaussian that doesn't mix spin up and spin down is broken #776, normalized so it satisfies W1 W1† + W2 W2† = I. This used to raise ValueError.
• test_bogoliubov_transform_spin_block_operator_algebra: checks U a_p† U⁻¹ = b_p† on the full unitary, for both a parity-preserving and a parity-flipping spin-up sector. This is the one that fails without the sign fix.
• test_bogoliubov_transform_spin_block_gaussian: eigenstate energies through the state-prep path.
• test_bogoliubov_transform_spin_mixing_gaussian_not_split: a transform that genuinely mixes spin should still go through the general path rather than getting split.

Outside the committed tests I ran a wider sweep (k = 2, 3, 4, real and complex, a range of seeds) cross-checking three things that each catch a sign error on their own: the operator algebra, whether U† H U is diagonal, and the eigenstate energies. Everything passes, including 27 of the parity-flipping cases. The full circuits/ suite passes and black/pylint are clean.

[google-cla]

Thanks for your pull request! It looks like this may be your first contribution to a Google open source project. Before we can look at your pull request, you'll need to sign a Contributor License Agreement (CLA).

View this failed invocation of the CLA check for more information.

For the most up to date status, view the checks section at the bottom of the pull request.

[gemini-code-assist]

Code Review

This pull request resolves issue #776 by introducing a new helper function _spin_blocks to correctly identify and split independent spin-sector transforms for both square and non-square (e.g., N x 2N) transformation matrices. This prevents non-square matrices from being misidentified as square spin-block-diagonal matrices, which previously led to slicing errors. Additionally, several unit tests and regression tests have been added to verify this behavior. There are no review comments, and I have no additional feedback to provide.

[blazingphoenix7]

@googlebot I signed it!

[mhucka]

Thank you for this work! Looks good to me. I requested Adam to review it too because this is a bit out of my depth.

[arettig]

Looks good and great catch on the parity flipping! Your tests look thorough as well.

I have one potential concern about the performance of the code - namely that the entire decomposition is run just to check the parity. I think there should be a way to do this more efficiently - I left a couple thoughts below. Let me know if it's not as straightforward as I'm thinking...

[arettig]

It would be nice to avoid running the entire decomposition just to check the parity if possible. I think there should be a straightforwardish way to check the parity of the transform just based on the transformation_matrix (maybe based on the singular values of the right_block or maybe the determinant?) that would speed this up a lot.