Hi @CatalinaAlbornoz -san,

Here are python codes to compare pennylane+lightning, pennylane+jax.jit, pennylane+jax.jit+jax.vmap, and qiskit-ml.

qiskit_machine_learning should be installed in advance.

The assumed use-case is QSVM with quantum kernel.

The number of circuits for kernel evaluation is 102x102 = 10^4.

[Main results]

pennylane+lightning : 26.5 sec.

pennylane+jax.jit : 12.4 sec.

pennylane+jax.jit+jax.vmap: 2.13 sec.

qiskit-ml : 1.34sec.

Sorry. My codes would not be beautiful.

(1) Pennylane+lightning

```
import pennylane as qml
from pennylane import numpy as np
from sklearn.svm import SVC
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
np.random.seed(1)
X, y = make_classification(n_features=2, n_redundant=0, n_informative=1, n_clusters_per_class=1, n_samples=1024)
# scaling the inputs is important since the embedding we use is periodic
scaler = StandardScaler().fit(X)
X_scaled = scaler.transform(X)
# scaling the labels to -1, 1 is important for the SVM and the
# definition of a hinge loss
y_scaled = 2 * (y - 0.5)
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, train_size=0.1)
n_qubits = 10
dev = qml.device("lightning.qubit", wires=n_qubits)
n_dim = len(X_train[0])
@qml.qnode(dev)
def kernel(x1, x2):
"""The quantum kernel."""
#S(x=x1)
for i in range(n_qubits):
qml.Hadamard(wires=[i])
qml.RZ(x1[i%n_dim], wires=[i])
#S^dagger(x=x2)
for i in range(n_qubits-1,-1,-1):
qml.RZ(-1*x2[i%n_dim], wires=[i])
qml.Hadamard(wires=[i])
return qml.probs(wires=range(n_qubits))
def kernel_matrix(A, B):
"""Compute the matrix whose entries are the kernel
evaluated on pairwise data from sets A and B."""
return np.array([[kernel(a, b)[0] for b in B] for a in A])
print(np.shape(X_train)[0],'x',np.shape(X_train)[0],'kernel_matrix','with', n_qubits,'qubits circuits')
%time kernel_matrix(X_train,X_train)
```

`102 x 102 kernel_matrix with 10 qubits circuits CPU times: user 26.4 s, sys: 78.1 ms, total: 26.5 s Wall time: 26.5 s`

(2) Pennylane+jax.jit

Note: jax might not work well on Windows.

```
# Added to silence some warnings.
from jax.config import config
config.update("jax_enable_x64", True)
import jax
from jax import device_put, jit
import pennylane as qml
from jax import numpy as np
from sklearn.svm import SVC
import numpy as vnp
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification
import matplotlib.pyplot as plt
vnp.random.seed(1)
X, y = make_classification(n_features=2, n_redundant=0, n_informative=1, n_clusters_per_class=1, n_samples=1024)
# scaling the inputs is important since the embedding we use is periodic
scaler = StandardScaler().fit(X)
X_scaled = scaler.transform(X)
# scaling the labels to -1, 1 is important for the SVM and the
# definition of a hinge loss
y_scaled = 2 * (y - 0.5)
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, train_size=0.1)
n_qubits = 10
dev = qml.device("default.qubit", wires=n_qubits)
n_dim = len(X_train[0])
X_train = device_put(X_train)
y_train = device_put(y_train)
@qml.qnode(dev, interface="jax")
def kernel(x1, x2):
"""The quantum kernel."""
#S(x=x1)
for i in range(n_qubits):
qml.Hadamard(wires=[i])
qml.RZ(x1[i%n_dim], wires=[i])
#S^dagger(x=x2)
for i in range(n_qubits-1,-1,-1):
qml.RZ(-1*x2[i%n_dim], wires=[i])
qml.Hadamard(wires=[i])
return qml.probs(wires=range(n_qubits))
jit_kernel = jax.jit(kernel)
def kernel_matrix(A, B):
"""Compute the matrix whose entries are the kernel
evaluated on pairwise data from sets A and B."""
return np.array([[jit_kernel(a, b).block_until_ready()[0] for b in B] for a in A])
print(np.shape(X_train)[0],'x',np.shape(X_train)[0],'kernel_matrix','with', n_qubits,'qubits circuits')
%time kernel_matrix(X_train,X_train)
```

`102 x 102 kernel_matrix with 10 qubits circuits CPU times: user 12 s, sys: 219 ms, total: 12.3 s Wall time: 12.4 s`

(3) Pennylane+jax.jit+jax.vmap

```
# Added to silence some warnings.
from jax.config import config
config.update("jax_enable_x64", True)
import jax
from jax import device_put, jit
import pennylane as qml
from jax import numpy as np
from sklearn.svm import SVC
import numpy as vnp
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification
import matplotlib.pyplot as plt
vnp.random.seed(1)
X, y = make_classification(n_features=2, n_redundant=0, n_informative=1, n_clusters_per_class=1, n_samples=1024)
# scaling the inputs is important since the embedding we use is periodic
scaler = StandardScaler().fit(X)
X_scaled = scaler.transform(X)
# scaling the labels to -1, 1 is important for the SVM and the
# definition of a hinge loss
y_scaled = 2 * (y - 0.5)
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, train_size=0.1)
n_qubits = 10
dev = qml.device("default.qubit", wires=n_qubits)
n_dim = len(X_train[0])
@qml.qnode(dev, interface="jax")
def kernel(x1, x2):
"""The quantum kernel."""
#S(x=x1)
for i in range(n_qubits):
qml.Hadamard(wires=[i])
qml.RZ(x1[i%n_dim], wires=[i])
#S^dagger(x=x2)
for i in range(n_qubits-1,-1,-1):
qml.RZ(-1*x2[i%n_dim], wires=[i]) #S(x)
qml.Hadamard(wires=[i])
return qml.probs(wires=range(n_qubits))
vectorized_kernel = jax.vmap(kernel)
jit_vectorized_kernel = jax.jit(vectorized_kernel)
# batching
result_0 = []
result_1 = []
for i in range(np.shape(X_train)[0]):
for k in range(np.shape(X_train)[0]):
result_0.append(X_train[i,:])
result_1.append(X_train[k,:])
x0_batch = np.array(result_0)
x1_batch = np.array(result_1)
print(np.shape(X_train)[0],'x',np.shape(X_train)[0],'kernel_matrix','with', n_qubits,'qubits circuits')
%time my_kernel_matrix = jit_vectorized_kernel(x0_batch,x1_batch).block_until_ready()
```

`102 x 102 kernel_matrix with 10 qubits circuits CPU times: user 4.47 s, sys: 1.14 s, total: 5.61 s Wall time: 2.13 s`

(4) qiskit_machine_learning

```
from sklearn.svm import SVC
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.datasets import make_classification
import matplotlib.pyplot as plt
X, y = make_classification(n_features=2, n_redundant=0, n_informative=1, n_clusters_per_class=1, n_samples=1024)
# scaling the inputs is important since the embedding we use is periodic
scaler = StandardScaler().fit(X)
X_scaled = scaler.transform(X)
# scaling the labels to -1, 1 is important for the SVM and the
# definition of a hinge loss
y_scaled = 2 * (y - 0.5)
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_scaled, train_size=0.1)
n_qubits = 10
n_dim = len(X_train[0])
from qiskit import BasicAer
from qiskit.circuit.library import ZZFeatureMap, ZFeatureMap, PauliFeatureMap
from qiskit.utils import QuantumInstance, algorithm_globals
from qiskit_machine_learning.algorithms import QSVC
from qiskit_machine_learning.kernels import QuantumKernel
n_reps = 1
input_data_dimension = n_dim
n_duplicate = int(n_qubits/input_data_dimension)
_X_train = np.tile(X_train, n_duplicate)
adhoc_feature_map = ZFeatureMap(feature_dimension=n_qubits, reps=1)
seed = 12345
adhoc_backend = QuantumInstance(BasicAer.get_backend('statevector_simulator'), shots=1, seed_simulator=seed, seed_transpiler=seed)
adhoc_kernel = QuantumKernel(feature_map=adhoc_feature_map, quantum_instance=adhoc_backend)
print(np.shape(X_train)[0],'x',np.shape(X_train)[0],'kernel_matrix','with', n_qubits,'qubits circuits')
%time adhoc_kernel.evaluate(_X_train,_X_train)
```

`102 x 102 kernel_matrix with 10 qubits circuits CPU times: user 672 ms, sys: 266 ms, total: 938 ms Wall time: 1.34 s`

I hope that these codes are helpful for pennylane lovers.