feat: support real machine gradient

.github/workflows/unittest.yml

@@ -22,7 +22,7 @@ jobs:
             run: python -m pip install torch torchvision torchaudio
 
         -   name: Install pyquafu
-            run: python -m pip install .
+            run: python setup.py develop
 
         -   name: Run unit tests
             run: pytest tests/

.gitignore

@@ -5,6 +5,7 @@
 *.so
 *.svg
 *temp.*
+temp/
 .DS_Store
 .idea
 .vscode

conftest.py

@@ -0,0 +1,17 @@
+import pytest
+
+
+def pytest_addoption(parser):
+    # Options for tests/quafu/algorithms/qnn_test.py
+    parser.addoption("--epoch", action="store", default=1, help="The number of epochs")
+    parser.addoption("--bsz", action="store", default=1, help="Batch size")
+
+
+@pytest.fixture
+def num_epochs(request):
+    return int(request.config.getoption("--epoch"))
+
+
+@pytest.fixture
+def batch_size(request):
+    return int(request.config.getoption("--bsz"))

quafu/algorithms/ansatz.py

@@ -145,16 +145,27 @@ class QuantumNeuralNetwork(Ansatz):
     """A Wrapper of quantum circuit as QNN"""
 
     # TODO(zhaoyilun): docs
-    def __init__(self, num_qubits: int, layers: List[Any], interface="torch"):
+    def __init__(
+        self, num_qubits: int, layers: List[Any], interface="torch", backend="sim"
+    ):
         """"""
         # Get transformer according to specified interface
         self._transformer = InterfaceProvider.get(interface)
         self._layers = layers
 
         # FIXME(zhaoyilun): don't use this default value
         self._weights = np.empty((1, 1))
+
+        self._backend = backend
         super().__init__(num_qubits)
 
+    def __call__(self, features):
+        """Compute outputs of QNN given input features"""
+        from .estimator import Estimator
+
+        estimator = Estimator(self, backend=self._backend)
+        return self._transformer.execute(self, features, estimator=estimator)
+
     def _build(self):
         """Essentially initialize weights using transformer"""
         for layer in self._layers:

quafu/algorithms/gradients/vjp.py

@@ -31,23 +31,36 @@ def _generate_expval_z(num_qubits: int):
 
 
 # TODO(zhaoyilun): support more measurement types
-def run_circ(circ: QuantumCircuit, params: Optional[List[float]] = None):
+# FIXME(zhaoyilun): remove backend
+def run_circ(
+    circ: QuantumCircuit,
+    params: Optional[List[float]] = None,
+    backend: str = "sim",
+    estimator: Optional[Estimator] = None,
+):
     """Execute a circuit
 
     Args:
         circ (QuantumCircuit): circ
         params (Optional[List[float]]): params
+        backend (str): backend
+        estimator (Optional[Estimator]): estimator
     """
     obs_list = _generate_expval_z(circ.num)
-    estimator = Estimator(circ)
+    if estimator is None:
+        estimator = Estimator(circ, backend=backend)
     if params is None:
         params = [g.paras for g in circ.parameterized_gates]
     output = [estimator.run(obs, params) for obs in obs_list]
     return np.array(output)
 
 
 # TODO(zhaoyilun): support more gradient methods
-def jacobian(circ: QuantumCircuit, params_input: np.ndarray):
+def jacobian(
+    circ: QuantumCircuit,
+    params_input: np.ndarray,
+    estimator: Optional[Estimator] = None,
+):
     """Calculate Jacobian matrix
 
     Args:
@@ -57,7 +70,8 @@ def jacobian(circ: QuantumCircuit, params_input: np.ndarray):
     batch_size, num_params = params_input.shape
     obs_list = _generate_expval_z(circ.num)
     num_outputs = len(obs_list)
-    estimator = Estimator(circ)
+    if estimator is None:
+        estimator = Estimator(circ)
     calc_grad = ParamShift(estimator)
     output = np.zeros((batch_size, num_outputs, num_params))
     for i in range(batch_size):
@@ -69,19 +83,38 @@ def jacobian(circ: QuantumCircuit, params_input: np.ndarray):
 
 
 def compute_vjp(jac: np.ndarray, dy: np.ndarray):
-    """compute vector-jacobian product
+    r"""compute vector-jacobian product
 
     Args:
         jac (np.ndarray): jac with shape (batch_size, num_outputs, num_params)
         dy (np.ndarray): dy with shape (batch_size, num_outputs)
+
+    Notes:
+        Suppose there are n inputs and m outputs in current node
+        Let x, y denote the inputs and outputs of current node, o denotes the final output
+        Essentially, jacobian is
+
+        .. math::
+            \begin{bmatrix}
+	    \frac{\partial y_1}{\partial x_1} & \cdots & \frac{\partial y_1}{x_n} \\
+	    \vdots & \ddots & \vdots \\
+	    \frac{\partial y_m}{\partial x_1} & \cdots & \frac{\partial y_m}{x_n}
+            \end{bmatrix}
+
+        `dy` is actually the vjp of dependent node
+
+        .. math:: \[ \frac{partial o}{partial y_1} \dots \frac{partial o}{partial y_m} \]
+
+        Therefore the vector jocobian product gets
+
+        .. math:: \[ \frac{partial o}{partial x_1} \dots \frac{partial o}{partial x_n} \]
     """
     batch_size, num_outputs, num_params = jac.shape
     assert dy.shape[0] == batch_size and dy.shape[1] == num_outputs
 
-    vjp = np.zeros((batch_size, num_params))
-
-    for i in range(batch_size):
-        vjp[i] = dy[i, :].T @ jac[i, :, :]
+    # Compute vector-Jacobian product using Einstein summation convention
+    #   the scripts simply mean 'jac-dims,dy-dims->vjp-dims'; so num_outputs is summed over
+    vjp = np.einsum("ijk,ij->ik", jac, dy)
 
     return vjp
 

quafu/algorithms/interface/torch.py

@@ -13,8 +13,11 @@
 # limitations under the License.
 """quafu PyTorch quantum layer"""
 
+from typing import Optional
+
 import numpy as np
 import torch
+from quafu.algorithms.estimator import Estimator
 
 from quafu import QuantumCircuit
 
@@ -28,14 +31,15 @@ def init_weights(shape):
         """Return torch gradient tensor with specified shape"""
         return torch.randn(*shape, requires_grad=True, dtype=torch.double)
 
-    # TODO(zhaoyilun): doc
+    # TODO(zhaoyilun): docstrings
     @staticmethod
     def execute(
         circ: QuantumCircuit,
         parameters: torch.Tensor,
         run_fn=run_circ,
         grad_fn=None,
         method="internal",
+        estimator: Optional[Estimator] = None,
     ):
         """execute.
 
@@ -45,11 +49,19 @@ def execute(
             grad_fn:
         """
 
-        kwargs = {"circ": circ, "run_fn": run_fn, "grad_fn": grad_fn}
+        kwargs = {
+            "circ": circ,
+            "run_fn": run_fn,
+            "grad_fn": grad_fn,
+            "estimator": estimator,
+        }
 
         if method == "external":
             return ExecuteCircuits.apply(parameters, kwargs)
         if method == "internal":
+            from ..ansatz import QuantumNeuralNetwork
+
+            assert isinstance(circ, QuantumNeuralNetwork)
             return ExecuteCircuits.apply(circ.weights, kwargs)
         raise NotImplementedError(f"Unsupported execution method: {method}")
 
@@ -61,11 +73,12 @@ class ExecuteCircuits(torch.autograd.Function):
     def forward(ctx, parameters, kwargs):
         ctx.run_fn = kwargs["run_fn"]
         ctx.circ = kwargs["circ"]
+        ctx.estimator = kwargs["estimator"]
         ctx.save_for_backward(parameters)
         parameters = parameters.numpy().tolist()
         outputs = []
         for para in parameters:
-            out = ctx.run_fn(ctx.circ, para)
+            out = ctx.run_fn(ctx.circ, para, estimator=ctx.estimator)
             outputs.append(out)
         outputs = np.stack(outputs)
         outputs = torch.from_numpy(outputs)
@@ -74,7 +87,7 @@ def forward(ctx, parameters, kwargs):
     @staticmethod
     def backward(ctx, grad_out):
         (parameters,) = ctx.saved_tensors
-        jac = jacobian(ctx.circ, parameters.numpy())
+        jac = jacobian(ctx.circ, parameters.numpy(), estimator=ctx.estimator)
         vjp = compute_vjp(jac, grad_out.numpy())
         vjp = torch.from_numpy(vjp)
         return vjp, None