Goal: Re-write documentation into Google style.

Reward: $100

The Symmer codebase is completely documented, but in an unofficial style. We would like the hacker to format our documentation into the Google style. The hacker should leave the descriptions at the beginning of each docstring as they are, but include necessary Arguments, Attributes, Returns etc. with relevant typing so that it is compliant with desired style. You may leave the description of arguments and returns blank, and the maintainers of Symmer will update this themselves. The winning PR will be awarded $100.

Description

QubitSubspaceManager.get_reduced_hamiltonian() has two arguments which default to None.

If no other value is given, the default causes an error.

Fix

Either

• Remove defaults
• Add logic to set a suitable value when None input

I'm running on a Mac M1, python 3.10 and no matter which version of Ray I install I always get the same error:

ModuleNotFoundError: No module named 'ray'

I am pretty sure that this is a more general ray issue which appears for some chips. Would it be possible to make ray option within symmer?

Add an automated subspace tool so that users do not need to understand any of the inner workings of symmer.

This problem lies at the heart of qubit tapering, where the assignment of eigenvalues to the identified symmetry generators places us in a so-called symmetry sector. The solution to the eigenvalue problem of interest (not limited to ground states, but excited states as well) will lie in one of these sectors, but it is far from trivial to determine which!

If a reference state can be identified which is consistent with the desired sector, then the problem becomes easy - simply measuring the symmetry generators in this state yields the necessary value assignments, thus placing us in the corresponding sector. Performing discrete optimization over the value assignments is hard - brute-forcing the solution scales exponentially with the size of the symmetry.

Currently only a single basis state (e.g. Hartree-Fock) is allowed as your input reference into stabilizer projection methods (tapering, CS-VQE etc.).

In github continuous integration check pip install poetry==1.4.0 to pip install poetry (when issues of poetry version 1.4.1 are fixed! Errors with debugpy)

Allow initialization from explicit representations, e.g. PauliwordOp.from_symplectic(.), PauliwordOp.from_dictionary(.), PauliwordOp.from_list(.), PauliwordOp.from_matrix(.) etc.

Currently Symmer has quite tight restrictions on python version:

[tool.poetry.dependencies]
python = ">=3.8 <3.11"

I'm going to suggest updating this to allow newer python versions to be included in the versioning (and then make it easier to use symmer in Nbed)

python = "^3.8"

Assess our performance against various quantum computing packages such as Qiskit, Cirq, Stim (the holy grail) etc.

Unitary partitioning is currently defined in our code only for M=2. However, in general we may have more cliques than this which we would like to be able to handle.

Could use m4ri library to speed up GF(2) calculations or potentially galois

This does not affect tapering, since the symmetry is of the full Hamiltonian, however in CS-VQE we identify a symmetry of a sub-Hamiltonian. Therefore, the basis of generators we select will influence the contextual terms that survive projection - we want to somehow maximise commutativity between the symmetry generators and the contextual Hamiltonian, since this will preserve maximal information of the problem in the contextual subspace.

e.g. suppose the symmetry generators are ZIII, IZII, IIZI and we have a contextual term XXXX. Then enforcing any of the stabilizers will result in the contextual term vanishing. Instead, if our generators are ZZII, IZZI, ZZZI then two out of the three preserve the contextual term, allowing us to drop two qubits-worth of quantum resource versus none.

Not clear what the optimal way of doing this is! So far, have experimented with defining Hamiltonian weighting functions that assign a score to each generator basis depending on the weight of the contextual terms it commutes with. Not limited to the pure Hamiltonian coefficients, for example have also tried a second-order response corrected approach based on https://arxiv.org/pdf/1406.4920.pdf.

• Simmer/Symmer
• Symplify

Minimal code to reproduce the bug:

Note the zero coefficient on the second rotation - this was causing one of the bugs in ContextualSubspace as the contextual operator was being mapped to zero.

Need to rotate onto +1 coefficient in the unitary partitioning step (seq_rot does this by default, hence not encountering the same issue). Can be achieved be conjugating with one of the other clique representatives (already tested this and is now giving the correct energies in the contextual subspace).

In the same python conda environment, line 28 of the below code (basically copied from one of your example notebooks) results in the file restarting multiple times and then eventually hitting an error when I use a python file rather than in a Jupyter notebook, where it does work. I am on a Mac if this matters. I have tried on two different computers and the problem persists.

The error goes down to symmer/operators/independent_op.py, line 238, in update_sector
with mp.Pool(mp.cpu_count()) as pool:

Which then raises an issue into the multiprocessing module, ultimately leading to a RuntimeError and then An attempt has been made to start a new process before the current process has finished its bootstrapping phase.

Code that reproduces the error below, the hamiltonian_data folder (copied from the repo) is in the same directory as the python/ipynb files

import numpy as np
import os
import json
from symmer.operators import PauliwordOp, QuantumState
from symmer.projection import QubitTapering

file_dir = os.getcwd()
ham_data_dir = os.path.join(file_dir, 'hamiltonian_data')

if not os.path.isdir(ham_data_dir):
raise ValueError('cannot find data dir')

filename = 'Be_STO-3G_SINGLET_JW.json'

if filename not in os.listdir(ham_data_dir):
raise ValueError('unknown file')

with open(os.path.join(ham_data_dir, filename), 'r') as infile:
data_dict = json.load(infile)
print("Successfully loaded data")

hf_state = QuantumState(np.asarray(data_dict['data']['hf_array'])) # Hartree-Fock state

hamiltonian_pauliword = PauliwordOp.from_dictionary(data_dict['hamiltonian'])

qubit_taper_object = QubitTapering(hamiltonian_pauliword)

H_taper = qubit_taper_object.taper_it(ref_state=hf_state)

print("completed")

Implement ADAPT-VQE in combination with CS-VQE. For example, using the contextual subspace greedy search heuristic of depth d (scales as O(N^{d+1})) this would work in the following way:

1. Define the ansatz pool (e.g. terms of the UCCSD(T) operator)
2. Choose the stabilizer pool (can just take the full set here, unless you have some additional intuition regarding the stabilizers and their importance)
3. Run over stabilizer relaxation combinations of size d
4. For each set of indices this will restrict the ansatz pool (i.e. perform the noncontextual projection over the pool)
5. Run qubit-ADAPT-VQE over the restricted ansatz pool
6. Repeat for each combination and choose the stabilizers that minimize the energy
7. Remove those stabilizers from the stabilizer pool and return to step 3
8. Repeat steps 3-8 until the stabilizer pool is exhausted
9. What remains will be the best stabilizer relaxation ordering and ansatze for the corresponding contextual subspaces

Description

Several places in Symmer use multiprocessing in the following way:

with mp.Pool(mp.cpu_count()) as pool: 

New process are created in one of two ways:

• Forking: default on linux and older MacOS systems
• Spawning: default on windows and newer MacOS
  Spawning a new process causes the top-level module to be re-imported to the new process which can cause infinite loops spawning new processes.

Minimal example

Running Symmer on MacOS Monterey (Apple M1 chip), I encountered this problem trying to run the code in notebook 2.2 as a local python script after cloning the Symmer repo to my local machine.

Suggested fixes

Ensuring all experiments run using Symmer are contained in the

if __name__ == "__main__":

idiom will prevent the top-level script being executed when it is imported to the spawned process stopping the infinite loop issue. This should then be documented in Symmer.

Alternatively, it's possible to force the new process to fork by using:

with mp.get_context("fork").Pool(mp.cpu_count()) as pool: 

This also fixes the issue without requiring the above if name == main idiom. However this may be bad practice (?). It would also need to be tested on Windows OS.

Goal: Write a faster to_sparse_matrix method for the PauliwordOp class in symmer.

Reward: $200

Details:

One of the core features of the symmer codebase is to provide a symbolic representation of the Pauli operators:

$$X =\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} $$

$$Y =\begin{pmatrix} 0 & -1i \\ 1i & 0 \end{pmatrix} $$

$$Z =\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} $$

$$I =\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} $$

The PauliwordOp class gives different functions that can be applied on these operators and tensor products of them. The individual operators can be written as:

from symmer import PauliwordOp

# dicitonary: {string: coefficient}
X = PauliwordOp.from_dictionary({"X":1})
Y = PauliwordOp.from_dictionary({"Y":1})
Z = PauliwordOp.from_dictionary({"Z":1})
I = PauliwordOp.from_dictionary({"I":1})

print(X)
print(Y)
print(Z)
print(I)

We can also represent tensor products of these Pauli operators as:

from symmer import PauliwordOp
# dicitonary: {string: coefficient}
XYZ = PauliwordOp.from_dictionary({"XYZ":1})
print(XYZ)

This object represents:

$$XYZ = X \otimes Y \otimes Z =\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \otimes \begin{pmatrix} 0 & -1i \\ 1i & 0 \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0& -1j& 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & +1j \\ 0 & 0 & 0 & 0 & +1j & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1j & 0 & 0\\ 0 & 0 & -1j & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & +1j & 0 & 0 & 0 & 0 \\ +1j & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -1j & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

and can be generated via:

from symmer import PauliwordOp
# dicitonary: {string: coefficient}
XYZ = PauliwordOp.from_dictionary({"XYZ":1})
print(XYZ.to_sparse_matrix.toarray())

Given the structure of Pauli operators $X, Y, Z, I$ rather than finding the kronecker product of each operator:

from scipy.linalg import kron
import numpy as np
from functools import reduce
X = np.array([[0, 1],
                       [1, 0]])
Y = np.array([[ 0., -1j],
                    [ 1j,  0.]])
Z = np.array([[1, 0],
                      [0, -1]])
I = np.eye(2)

XYZ = reduce(kron, [X,Y,Z])
print(XYZ)

We use a to_sparse_matrix method, that determines the indices of the output matrix from the structure of the operators. This increases the speed in finding the matrix respresentation of the operator over manually taking the individual kronecker products. However, there are still improvements that can be made hence this pull request. We want someone to speed up this functionality.

One approach is outlined in a qiskit issue: Qiskit/qiskit#8772 , that points to: https://github.com/chetmurthy/qrusty. This codebase is written in rust. It would be great to have someone convert this rust implementation into a python version for symmer.

Alternatively, if someone can find a different method that is faster than the current one this would also be a valid solution!

The winning PR will be awarded $200.

The current CS-VQE implementation prioritises diagonal Hamiltonian terms (hence why the noncontextual energy aligns closely with the Hartree-Fock energy). This is not necessary, for example we might maximise anticommuting contributions - the noncontextual energy will be (considerably) worse in this case, but how does this affect the resulting quantum corrections?

Goal: Optimize MPO construction in Symmer.

Reward: $200

A subroutine of qubit subspace techniques is the identification of the correct symmetry sector - an assignment of $\pm1$ eigenvalues to the chosen stabilizers over which we project to yield a reduced subspace. Naively this problem scales exponentially in the number of stabilizers one wishes to enforce; however, if you have access to a reference state with non-zero overlap with the ground state, the eigenvalue-assignment problem becomes trivial - 'measure' each symmetry generator in the reference state and we get the assignment for free. This poses a problem - how do we identify such a reference state? For quantum chemistry application we may choose the Hartree-Fock state, but in general it is not straightforward.

One approach is to use the Density Matrix Renormalization Group (DMRG), that allows you to construct approximate ground states, whose overlap is dependent on the bond dimension of an associated Matrix Product Operator (MPO).

This is implemented in symmer.approximate, however there are optimizations that can be made to speed-up the construction of MPOs, for example it could be parallelised.

That is the goal of this bounty.

The success of CS-VQE is highly sensitive to the set of stabilizers one wishes to enforce.

Stabilizers should be selected such that they maximize commutation with respect to some weighting function so that optimal information is preserved under projection onto the contextual subspace.