Raise error for invalid `section_size` in `synth_cnot_count_full_pmh`

[Tarun-Kumar07]

Summary

This PR enforces the validation that section_size must be less than or equal to num_qubits for proper functioning of the synth_cnot_count_full_pmh algorithm, as described in Efficient Synthesis of Linear Reversible Circuits. It fixes #12106

Description of change

• Added validation to enforce section_size <= number of qubits.
• Defaulted section_size to proportional to np.log(num_qubits).
  • This leads to decrease in depth of synthesised circuits compared to previous default value of 2

[qiskit-bot]

Thank you for opening a new pull request.

Before your PR can be merged it will first need to pass continuous integration tests and be reviewed. Sometimes the review process can be slow, so please be patient.

While you're waiting, please feel free to review other open PRs. While only a subset of people are authorized to approve pull requests for merging, everyone is encouraged to review open pull requests. Doing reviews helps reduce the burden on the core team and helps make the project's code better for everyone.

One or more of the the following people are requested to review this:

• @Qiskit/terra-core

[CLAassistant]

All committers have signed the CLA.

[qiskit-bot]

Thank you for opening a new pull request.

Before your PR can be merged it will first need to pass continuous integration tests and be reviewed. Sometimes the review process can be slow, so please be patient.

While you're waiting, please feel free to review other open PRs. While only a subset of people are authorized to approve pull requests for merging, everyone is encouraged to review open pull requests. Doing reviews helps reduce the burden on the core team and helps make the project's code better for everyone.

One or more of the the following people are requested to review this:

• @Qiskit/terra-core

[coveralls]

Pull Request Test Coverage Report for Build 8785363931

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

• For more information on this, see Tracking coverage changes with pull request builds.
• To avoid this issue with future PRs, see these Recommended CI Configurations.
• For a quick fix, rebase this PR at GitHub. Your next report should be accurate.

Details

• 3 of 3 (100.0%) changed or added relevant lines in 1 file are covered.
• 1 unchanged line in 1 file lost coverage.
• Overall coverage increased (+0.03%) to 89.295%

---

Files with Coverage Reduction	New Missed Lines	%
crates/qasm2/src/lex.rs	1	92.62%
<!--	Total:	1

Totals
Change from base Build 8783263005:	0.03%
Covered Lines:	60369
Relevant Lines:	67606

---

💛 - Coveralls

[Tarun-Kumar07]

Thank you, @ShellyGarion. I'll experiment with the default value. Do we have a documentation repository where we document such experiments for future reference?

[ShellyGarion]

Thank you, @ShellyGarion. I'll experiment with the default value. Do we have a documentation repository where we document such experiments for future reference?

Please summarize your results in the PR comments.

[Tarun-Kumar07]

Hey @ShellyGarion,

I'm not sure why the code below is failing when I provide different section_size values. The unitary matrices of the circuits don't match.

import numpy.testing as np
import qiskit.quantum_info as qi

from qiskit.synthesis import synth_cnot_count_full_pmh
from qiskit.synthesis.linear import random_invertible_binary_matrix

mat = random_invertible_binary_matrix(6, seed=1234)
qc_4 = synth_cnot_count_full_pmh(mat, section_size=4)
op_4 = qi.Operator(qc_4)

qc_2 = synth_cnot_count_full_pmh(mat, section_size=2)
op_2 = qi.Operator(qc_2)

np.assert_equal(op_4.data, op_2.data) # Fails

For the given binary matrix of shape (6,6), something seems off with section_size = 4. I repeated the above code for different section_size values and compared their unitary matrices, the results are below:

	1	2	3	4	5	6
1	:white_check_mark:	:white_check_mark:	:white_check_mark:	:x:	:white_check_mark:	:white_check_mark:
2	:white_check_mark:	:white_check_mark:	:white_check_mark:	:x:	:white_check_mark:	:white_check_mark:
3	:white_check_mark:	:white_check_mark:	:white_check_mark:	:x:	:white_check_mark:	:white_check_mark:
4	:x:	:x:	:x:	:white_check_mark:	:x:	:x:
5	:white_check_mark:	:white_check_mark:	:white_check_mark:	:x:	:white_check_mark:	:white_check_mark:
6	:white_check_mark:	:white_check_mark:	:white_check_mark:	:x:	:white_check_mark:	:white_check_mark:

Do you have any insight into why this is happening?

Thanks!

[alexanderivrii]

@Tarun-Kumar07, thanks for spotting this bug - now the issue is significantly more problematic than what I originally thought. A quick experiment on my end (varying the number of qubits and considering 100 random invertible binary matrices for each) gives the following output:

For nq = 2: bad_section_sizes = set()
For nq = 3: bad_section_sizes = set()
For nq = 4: bad_section_sizes = set()
For nq = 5: bad_section_sizes = {3}
For nq = 6: bad_section_sizes = {4}
For nq = 7: bad_section_sizes = {4, 5}
For nq = 8: bad_section_sizes = {3, 5, 6}
For nq = 9: bad_section_sizes = {5, 6, 7}

Would you be interested to take a deeper look both at the paper and the code and see if you can fix the problem?

[alexanderivrii]

For reference, here is my python script (based on yours):

from qiskit.quantum_info import Operator
from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import LinearFunction
from qiskit.synthesis import synth_cnot_count_full_pmh
from qiskit.synthesis.linear import random_invertible_binary_matrix

for nq in range(2, 10):
    bad_section_sizes = set()

    for seed in range(100, 200):
        mat = random_invertible_binary_matrix(nq, seed=seed)

        expected_qc = QuantumCircuit(nq)
        expected_qc.append(LinearFunction(mat), range(nq))
        expected_op = Operator(expected_qc)

        for section_size in range(1, nq+1):
            qc = synth_cnot_count_full_pmh(mat, section_size=section_size)
            op = Operator(qc)
            ok = op.equiv(expected_op)
            # print(f"section size is {section_size}: {ok}")
            if not ok:
                bad_section_sizes.add(section_size)

    print(f"For {nq = }: {bad_section_sizes = }")

[Tarun-Kumar07]

@alexanderivrii , I will investigate further and let you know

[ShellyGarion]

In a previous version of this code there was the following comment on section_size: section_size must be a factor of the number of qubits. (it was removed in #12107). However, it seems that it's not always necessary, since section_size = 2 works also when the number of qubits is odd. So we are wondering what should be the conditions on section_size. And of course, if an invalid value of section_size is provided, the code should raise an error, and not produce a wrong circuit.

[ShellyGarion]

@Tarun-Kumar07 - are you still working on this ?

[Tarun-Kumar07]

Hi @ShellyGarion

I'm currently tied up with some personal matters and will be able to resume work after June 15th. Is that okay with you?

Thanks for understanding!

[Abdalla01001]

@ShellyGarion @alexanderivrii The issue in the code stems from how the algorithm handles the section_size parameter. The algorithm is designed to divide the input matrix into sections and process these sections iteratively. However, when the section_size is larger than the matrix size, there are no sections to iterate over, causing the algorithm to fail. I think we can try to solve this by virtual padding, If the section_size exceeds the matrix size (n), we can virtually pad the matrix with additional columns of zeros to make it divisible by the section_size. This padding is "virtual" in the sense that it's not explicitly added to the matrix data structure, but rather accounted for in the algorithm's logic.

I want you hopefully to reproduce this edited version of function synth_cnot_count_full_pmh and review Its reliability :

import numpy as np
from qiskit.circuit import QuantumCircuit
from qiskit.exceptions import QiskitError

def synth_cnot_count_full_pmh(
    state: list[list[bool]] | np.ndarray[bool], section_size: int = 2
) -> QuantumCircuit:
    """Synthesize linear reversible circuits using the Patel-Markov-Hayes algorithm with virtual padding."""

    if not isinstance(state, (list, np.ndarray)):
        raise QiskitError(
            "state should be of type list or numpy.ndarray, "
            "but was of the type {}".format(type(state))
        )
    state = np.array(state)
    num_qubits = state.shape[0]

    # Virtual Padding
    if section_size > num_qubits:
        padding_size = section_size - (num_qubits % section_size)
        state = np.hstack((state, np.zeros((num_qubits, padding_size), dtype=bool)))

    # Synthesize lower triangular part
    [state, circuit_l] = _lwr_cnot_synth(state, section_size, num_qubits)
    state = np.transpose(state)

    # Synthesize upper triangular part
    [state, circuit_u] = _lwr_cnot_synth(state, section_size, num_qubits)
    circuit_l.reverse()
    for i in circuit_u:
        i.reverse()

    # Convert the list into a circuit of C-NOT gates (removing virtual padding)
    circ = QuantumCircuit(num_qubits)  # Circuit size is the original num_qubits
    for i in circuit_u + circuit_l:
        # Only add gates if both control and target are within the original matrix
        if i[0] < num_qubits and i[1] < num_qubits:
            circ.cx(i[0], i[1])

    return circ

def _lwr_cnot_synth(state, section_size, num_qubits):
    """Helper function for the Patel-Markov-Hayes algorithm with virtual padding."""

    circuit = []
    cutoff = 1

    # Iterate over column sections (including padded sections)
    for sec in range(1, int(np.ceil(state.shape[1] / section_size)) + 1):
        # Remove duplicate sub-rows in section sec
        patt = {}
        for row in range((sec - 1) * section_size, num_qubits):
            sub_row_patt = copy.deepcopy(
                state[row, (sec - 1) * section_size : sec * section_size]
            )
            if np.sum(sub_row_patt) == 0:
                continue
            if str(sub_row_patt) not in patt:
                patt[str(sub_row_patt)] = row
            else:
                state[row, :] ^= state[patt[str(sub_row_patt)], :]
                circuit.append([patt[str(sub_row_patt)], row])

        # Use gaussian elimination for remaining entries in column section
        # Modified loop range
        for col in range((sec - 1) * section_size, min(sec * section_size, num_qubits)):
            # Check if 1 on diagonal
            diag_one = 1
            if state[col, col] == 0:
                diag_one = 0
            # Remove ones in rows below column col
            for row in range(col + 1, num_qubits):
                if state[row, col] == 1:
                    if diag_one == 0:
                        state[col, :] ^= state[row, :]
                        circuit.append([row, col])
                        diag_one = 1
                    state[row, :] ^= state[col, :]
                    circuit.append([col, row])
                # Back reduce the pivot row using the current row
                if sum(state[col, :] & state[row, :]) > cutoff:
                    state[col, :] ^= state[row, :]
                    circuit.append([row, col])

    return [state, circuit]

• Virtual Padding: If the section_size is larger than the number of qubits (num_qubits), the matrix is padded with zeros to the right using np.hstack.

• Modified Loop: The loop in _lwr_cnot_synth now iterates up to state.shape[1] (the new padded size) to cover all sections.

• Handling Zero Sub-rows: The code now checks if a sub-row is all zeros (np.sum(sub_row_patt) == 0) and skips it if so.

• Removing Padding: When creating the circuit, we only add CNOT gates if both the control and target qubits are within the original num_qubits, effectively discarding operations on the padded columns.

• Modified Loop Range: The range of the inner loop in _lwr_cnot_synth is modified to range((sec - 1) * section_size, min(sec * section_size, num_qubits)). This ensures that the loop only iterates over the original columns, even if the section_size is larger than the matrix size.

[ShellyGarion]

@Abdalla01001 - thanks for your suggestion. Could you perhaps verify the correctness of your suggested code? See the following comment of @Tarun-Kumar07 https://github.com/Qiskit/qiskit/pull/12166#issuecomment-2088040391

[Abdalla01001]

@ShellyGarion Yeah: For nq = 2: bad_section_sizes = set() For nq = 3: bad_section_sizes = set() For nq = 4: bad_section_sizes = set() For nq = 5: bad_section_sizes = set() For nq = 6: bad_section_sizes = set() For nq = 7: bad_section_sizes = set() For nq = 8: bad_section_sizes = set() For nq = 9: bad_section_sizes = set()

[ShellyGarion]

thanks @Abdalla01001 for checking this. Could you also check that with your code, for large number of qubits (say, 20, 50, 100, 200, 500, 1000), the value of section_size = np.log(num_qubits) is better than the current default value (=2) ?

[Abdalla01001]

@ShellyGarion It works correctly, but I avoid waiting for 500 and 1000 qubits because of my local resources, It can take above 2 hours, but you can check it! I made a comparison with various number of qubits for the two sec_size default values, It says that the log-based will be less efficient because the CX count in the most cases is higher:

num_qubits_list = [20, 50, 100, 200]
results = {}

for num_qubits in num_qubits_list:
    cnot_counts_default = []
    cnot_counts_log = []

    for _ in range(10):  # Repeat 10 times for each num_qubits
        mat = random_invertible_binary_matrix(num_qubits)

        # Default
        cnot_gates_default = synth_cnot_count_full_pmh(mat, section_size=2)
        cnot_counts_default.append(len(cnot_gates_default))

        # Log-based
        section_size_log = round(np.log2(num_qubits))
        cnot_gates_log = synth_cnot_count_full_pmh(mat, section_size=section_size_log)
        cnot_counts_log.append(len(cnot_gates_log))

    results[num_qubits] = {
        "default_avg": np.mean(cnot_counts_default),
        "log_avg": np.mean(cnot_counts_log),
    }

print(results)

results:

{20: {'default_avg': 294.7,   'log_avg': 297.6},
50: {'default_avg': 2140.1,   'log_avg': 2227.6}, 
100: {'default_avg': 8726.5,   'log_avg': 9109.0},
200: {'default_avg': 35049.0,   'log_avg': 36576.6}}

[Tarun-Kumar07]

Hi @Abdalla01001, Could you please take over the PR and this issue ? Let me know if you need some help. Thanks !!

[Abdalla01001]

@Tarun-Kumar07 sure!

[Abdalla01001]

@ShellyGarion @alexanderivrii Any updates?

[ShellyGarion]

@ShellyGarion @alexanderivrii Any updates?

@Abdalla01001 - I've already assigned you in the issue #12106, so feel free to open a new PR, and then we'll close this one. You can add there your suggested code, as well as tests to validate the correctness of the circuit synthesis. Thanks!

[Abdalla01001]

@ShellyGarion Done! #12712

[ShellyGarion]

As discussed above, I'm closing this PR in favor of #12712

[Cryoris]

In the paper they provide an upper bound on the number of row ops (-> number of CX ops)

If we plot the optimal section size $m$ and fit it with a model $m = \alpha\log_2(n)$ we can try to fit the parameter. From that it seems that:

1. up to 100 qubits section_size=2 is best
2. afterwards $\alpha\approx 0.56$ should perform well

In the above plot suggestion is the function max(2, alpha * np.log2(n)) which we could use as default value.