apply_gate() results mismatch with qiskit?

[wcqc]

I initialized the state in cuquantum (2nd code block below) with the state from Qiskit after applying the first HGate, then applied the CXGate to the state in cuquantum, somehow the final resulting state differs from the one obtained from Qiskit? Is there anything not looking right in the 2nd code block blow, anything obviously incorrect/missing etc?

expected: tensor([[0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.0000+0.j,
         0.0000+0.j]], device='cuda:0'), result: tensor([[0.7071+0.j, 0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.0000+0.j, 0.0000+0.j, 0.0000+0.j,
         0.0000+0.j]], device='cuda:0')

from qiskit import QuantumCircuit
import qiskit.quantum_info as qi
from qiskit.extensions import CXGate

cx = CXGate()
print(cx.to_matrix())

qc = QuantumCircuit(3)
state = qi.Statevector.from_instruction(qc)
print(state)
qc.h(0)
state = qi.Statevector.from_instruction(qc)
print(state)
qc.cx(0, 1)
state = qi.Statevector.from_instruction(qc)
print(state)

import torch
from cuquantum import custatevec as cusv
from cuquantum import cudaDataType, ComputeType

C_DTYPE = torch.complex64

import torch

nSVs = 1
nIndexBits = 3
svSize = (1 << nIndexBits)
svStride = svSize
adjoint = 0

targets = [1]
nTargets = len(targets)
controls = [0]
nControls = len(controls)

matrixIndices = [0]
nMatrices = len(matrixIndices)

print('bsz = 1')
d_svs = torch.tensor([[0.70710678 + 0.j, 0.70710678 + 0.j, 0. + 0.j,
                       0. + 0.j, 0. + 0.j, 0. + 0.j,
                       0. + 0.j, 0. + 0.j]], dtype=C_DTYPE).cuda()
d_svs_res = torch.tensor([[0.70710678 + 0.j, 0. + 0.j, 0. + 0.j,
                           0.70710678 + 0.j, 0. + 0.j, 0. + 0.j,
                           0. + 0.j, 0. + 0.j],
                          ], dtype=C_DTYPE).cuda()
d_matrices = torch.tensor([[1. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
                           [0. + 0.j, 0. + 0.j, 0. + 0.j, 1. + 0.j],
                           [0. + 0.j, 0. + 0.j, 1. + 0.j, 0. + 0.j],
                           [0. + 0.j, 1. + 0.j, 0. + 0.j, 0. + 0.j]], dtype=C_DTYPE).cuda()

###################################################################################a

# cuStateVec handle initialization
handle = cusv.create()
# check the size of external workspace
extraWorkspaceSizeInBytes = cusv.apply_matrix_batched_get_workspace_size(
    handle, cudaDataType.CUDA_C_32F, nIndexBits, nSVs, svStride,
    cusv.MatrixMapType.BROADCAST, matrixIndices, d_matrices.data_ptr(),
    cudaDataType.CUDA_C_32F, cusv.MatrixLayout.ROW, adjoint, nMatrices,
    nTargets, nControls,
    ComputeType.COMPUTE_32F)

workspace_ptr = 0

# apply gate
cusv.apply_matrix_batched(
    handle, d_svs.data_ptr(), cudaDataType.CUDA_C_32F, nIndexBits, nSVs, svStride,
    cusv.MatrixMapType.BROADCAST, matrixIndices, d_matrices.data_ptr(),
    cudaDataType.CUDA_C_32F, cusv.MatrixLayout.ROW, adjoint, nMatrices,
    targets, nTargets, controls, 0, nControls,
    ComputeType.COMPUTE_32F, workspace_ptr, extraWorkspaceSizeInBytes)

# destroy handle
cusv.destroy(handle)

# check result
if not torch.allclose(d_svs_res, d_svs):
    print(f'expected: {d_svs_res}, result: {d_svs}')
    raise ValueError("results mismatch")
print(d_svs, d_svs_res)
print("test passed")

[mtjrider]

Hi @wcqc

I've annotated your code with ↑↓ to highlight the problematic portions of your example.
You haven't specified the 2-qubit gate contiguously (concatenated by row). Based on your current definition, each row of the 2-qubit gate matrix is being treated as an individual gate. You are targeting qubit [1] with control qubit [0] using a gate defined by the first row of the matrix [1. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j] => [[1. + 0.j, 0. + 0.j], [0. + 0.j, 0. + 0.j]].
 

...
↓↓↓↓↓↓↓↓↓↓↓↓↓↓
targets = [1]
↑↑↑↑↑↑↑↑↑↑↑↑↑↑
nTargets = len(targets)
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
controls = [0]
↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑
nControls = len(controls)

matrixIndices = [0]
nMatrices = len(matrixIndices)

print('bsz = 1')
d_svs = torch.tensor([[0.70710678 + 0.j, 0.70710678 + 0.j, 0. + 0.j,
                       0. + 0.j, 0. + 0.j, 0. + 0.j,
                       0. + 0.j, 0. + 0.j]], dtype=C_DTYPE).cuda()
d_svs_res = torch.tensor([[0.70710678 + 0.j, 0. + 0.j, 0. + 0.j,
                           0.70710678 + 0.j, 0. + 0.j, 0. + 0.j,
                           0. + 0.j, 0. + 0.j],
                          ], dtype=C_DTYPE).cuda()
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
d_matrices = torch.tensor([[1. + 0.j, 0. + 0.j, 0. + 0.j, 0. + 0.j],
                           [0. + 0.j, 0. + 0.j, 0. + 0.j, 1. + 0.j],
                           [0. + 0.j, 0. + 0.j, 1. + 0.j, 0. + 0.j],
                           [0. + 0.j, 1. + 0.j, 0. + 0.j, 0. + 0.j]], dtype=C_DTYPE).cuda()
↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑
...

To fix this, you have a couple of options:

1. treat the 2-qubit gate as a dense matrix, making it contiguous by concatenating rows. Change targets so that it points to the control-ordered qubits, and switch controls = 0; nControls = 0.
2. implement the CX operation natively with the batched apply matrix API.

To illustrate (2), I've written a gist and posted it here. I wrote it somewhat hastily, so let me know if you run into problems.

Edit: the Python gist should print this:

calculated:
tensor([[0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.0000+0.j,
         0.0000+0.j]], device='cuda:0')
expected:
tensor([[0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.7071+0.j, 0.0000+0.j, 0.0000+0.j, 0.0000+0.j,
         0.0000+0.j]], device='cuda:0')
test passed

[wcqc]

Hi @mtjrider,

I've tried both of the 2 options and now the test passed, and just to confirm: when the gate is not decomposed (like in option 1), we always have to encode both "controls" and "targets" in targets instead of separating them into targets and controls (like in option 2).

Is option 2 preferred over 1 in terms of efficiency?

[mtjrider]

@wcqc yes. When the gate is not decomposed, you must provide both controls and targets to targets=[] where you’ve taken care in the relative qubit index ordering to preserve the role of the controls.

We tend to prefer option 2 as it provides more composability— gates can be fused, reordered, etc..

The assumption we make in the apply_matrix API is that the gate is numerically dense. We have other APIs for when the gate takes a special form (e.g.) Pauli matrices, diagonal gates, rotations, etc.

I’ll add that there are cases when fusing the control in this fashion can be advantageous. It allows you to fuse 1-qubit gates on either qubit before and after the CX. Fusion is useful because it increases arithmetic intensity and reduces the number of API calls required by a gate sequence.

[wcqc]

One additional question regarding gate fusion: when the gates are not decomposed into single-qubit gates, gate fusion can still kick in (when the cuquantum API is called) as long as there are fusible gates in the circuit? When you say " there are cases when fusing the control in this fashion can be advantageous", this is pointing to cases where there are extra gains on top of the fusions applicable in a non-decomposed circuit (e.g., a circuit mixed with 1- and 2-qubit gates)? Or do you actually mean to take advantage of any cuquantum gate fusion, all gates have to be decomposed into single-qubit gates?

Edit: In order to use multi-GPU simulation with cuStateVec, is it always necessary to decompose all gates into single-qubit ones?

[mtjrider]

@wcqc I see that I may have confused you.

Gate fusion is not done automatically by these APIs. You must implement this yourself as we don't currently provide a transpiler. There are many other open-source solutions that address this problem. The only "exception" is the cuQuantum Appliance container, which is a complete quantum circuit simulation framework.

Very loosely speaking, I meant this-- given two gate operations (A and B) on potentially different qubit sets, assume the following:

1. The execution time of A is $t_A$ and the execution time of B is $t_B$.
2. The overhead for calling any gate is constant and takes time $T_o$.

If we execute sequentially, the total runtime is $t_r = 2 T_o + t_A + t_B$. Depending on the hardware running the gate application, fusing the gates will not increase the runtime by relying on the physical hardware to accommodate the increased computational intensity. In the best case, the new runtime would be: $t_r = T_o + t_C$, where $t_C$ is the runtime for the fused operation.

We can avoid API call overhead, and potentially shorten the effective runtime (the difference being $t_d = t_A + t_B - t_C$ ). This is only possible if the hardware running the operation can accommodate the larger, fused operation.

A final note: (2) is not always true but simplifies my description.

[wcqc]

@mtjrider Somehow I got confused by reading the gate fusion section in the doc (https://docs.nvidia.com/cuda/cuquantum/23.06.1/custatevec/overview.html) and believed the API does gate fusion. :)

[mtjrider]

Ah, yes. That document describes how you'd implement gate fusion using our apply_matrix API.
Is there any specific language that lead to the confusion?

[mtjrider]

Is this a decomposition/fusion of the bell state circuit?

@wcqc let me clarify something which is alluded to by our multi-GPU documentation.

$$ \begin{align} \mathrm{C} X &= \ket{0} \bra{0} \otimes I + \ket{1} \bra{1} \otimes X \\ \dots &= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix} \\ \Rightarrow \mathrm{C}X \ket{\psi} &= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \psi_{00} \\ \psi_{01} \\ \psi_{10} \\ \psi_{11} \end{bmatrix} = \begin{bmatrix} \psi_{00} \\ \psi_{01} \\ \psi_{11} \\ \psi_{10} \end{bmatrix} \end{align} $$

Please note that $\ket{\psi}$ represents a 2-qubit state, and I've written the indices of the state vector in binary (e.g.) $\psi_{i}$ with $i \in \left[ 2^{n {\rightarrow 2}} \right]$ $\rightarrow$ $\psi_{q_{0},q_{1}}$ with $q_{i} \in \left[ 2 \right]$. (1)

With these observations, $\mathrm{C}X \ket{\psi} \equiv X \psi_{q_{0} {\rightarrow 1},\ q_{1}}$. This doesn't represent a decomposition, per se. It's one prescription for handling control sequences with arbitrary gates without forming the entire 2-qubit gate. In this way, we reduce arithmetic intensity for whatever penalty is implied by the index slicing.

This means that you can implement a 2-qubit Bell state via:

import numpy as np


dtype = np.complex64
val = np.sqrt(2)/2
H = np.array(
    [[val, val],
     [val, -val]],
    dtype=dtype
)
X = np.array(
    [[0, 1],
     [1, 0]],
    dtype=dtype
)

# use index slicing
psi = np.zeros(4, dtype=dtype).reshape(2,2)
psi[0, 0] = 1
## apply H to the most significant qubit in the array
for qi in range(2):  # the local dimension of a single site
    psi[:, qi] = H @ psi[:, qi]  # this operation targets qubit 0
## apply X on qubit 1 controlled by qubit 0
psi[1, :] = X @ psi[1, :]
## reshape the state to a vector
psi = psi.reshape(4)


# form 2-qubit operators directly
phi = np.zeros(4, dtype=dtype)
phi[0] = 1
H_2qubit = np.kron(H, np.eye(2))
CX = (
    np.kron(np.outer(*[[1, 0]]*2), np.eye(2)) +
    np.kron(np.outer(*[[0, 1]]*2), X)
)
phi = CX @ H_2qubit @ phi

if np.allclose(psi, phi):
    print("these are equal")
else:
    print("these are not equal")

(1) $[n] \coloneqq [0,1,\dots,n-1]$

[wcqc]

The above CX matrix assumes big-endian, and cuquantum uses little-endian and in the latter case the above "decomposition" would not apply assuming little-endianness (because the right bottom corner matrix would not be pauli-X)? Would the above best be done by a transpiler or is it just a basic and general linear algebra op that can be applied to handling control sequences with arbitrary gates without forming the entire 2-qubit gate, and should be done manually for best efficiency when using cuquantum?

[mtjrider]

I'm somewhat confused by your question. When a control-qubit is involved, there is no operational decomposition.
The effect of a control is to introduce a (potentially complex) index slicing scheme on the state vector based on the product space of local sites (qubits).

In short, it does not matter if a state vector is expressed in big-/little- endian when considering something like C-X. The only thing that matters-- which qubit index is the control, what is the "value" of that control, and what qubit index is the target.

More succinctly: $\mathrm{C}X \ket{\psi} \equiv X \psi_{q_{0} {\rightarrow 1},\ q_{1}}$, can be generalized to $\mathrm{C}X \ket{\psi} \equiv X \psi_{ q_{0},\ q_{1},\ \cdots\ \ q_{i} \rightarrow c,\ \cdots\ q_{n-1}}$ where $c$ is the control-value ($0$, $1$ for qubits ...) ... without regard to endianness (the choice of bit-wise significance in $q_{i}$ is arbitrary).

The end-result is that you need only supply our API with the dense operator (in this example, $X$), and the control qubit index + value. If multiple controls/values are involved, a list may be provided.

If this does not make sense, can you provide further clarification regarding the confusion?

[wcqc]

In the very top example when this thread first started, where the 2-qubit dense CXgate is specified as the matrix/gate to be applied to the 3-qubit statevec, cuQuantum would automatically "expand" the 2-qubit dense CXgate matrix (4x4) to obtain the dense CXgate (8x8) for 3-qubits under the hood, and this expansion would happen automatically for any call of the api involving say a multi-qubit dense gate matrix supplied as long as we do the correct targets and controls spec (by putting everything in targets)?

In the second case where the bell state is obtained by an Hgate followed by an Xgate, the expansion would be only on the dense X; the api knows this as X is a single-qubit gate?

And the API distinguishes the above two cases by looking at how the targets and controls are specified, eg. in the first case both are in targets?

[mtjrider]

Hi @wcqc thanks for your clarification. I apologize if my previous responses caused confusion.

Let me simplify things a bit. There are three cases we've looked at:

Note: a "dense m-qubit gate" is a non-sparse matrix with dimensions $2^{m} \times 2^{m}$

1. Passing $\mathrm{C}X$ as a dense 2-qubit gate to apply_matrix.
2. Embedding $\mathrm{C}X$ into a dense 3-qubit gate, and passing it to apply_matrix.
3. Passing only $X$ to apply_matrix with appropriately defined targets and controls.

Item (3) is the simplest and most important case to understand. If you know that the gate operation is fundamentally representable with a dense 1-qubit gate, controls, and targets, then you do not need to embed that 1-qubit gate into a larger m-qubit gate unless you are implementing gate-fusion or some other optimization.

apply_matrix -- or any other API from custatevec -- currently does not implement gate-fusion, etc., for you. You must implement it yourself. One reason to embed $\mathrm{C}X$ in a larger gate is to then fuse it with other gates.

When you pass a dense 2-qubit gate to apply_matrix and specify controls, it is interpreted as follows:

... meaning, $a \equiv \mathrm{C}X$, and this 2-qubit gate is controlled by a third qubit.