diff --git a/.github/workflows/deploy_ghpages.yml b/.github/workflows/deploy_ghpages.yml index b61cf7da..70bce3b0 100644 --- a/.github/workflows/deploy_ghpages.yml +++ b/.github/workflows/deploy_ghpages.yml @@ -22,6 +22,12 @@ jobs: with: path: ~/.cache/pip key: deploy_ghpages.yml + - name: Copy feature maps + run: | + cp -a ./symb_statevectors/. ./examples/ERP/symb_statevectors + cp -a ./symb_statevectors/. ./examples/MI/symb_statevectors + cp -a ./symb_statevectors/. ./examples/other_datasets/symb_statevectors + cp -a ./symb_statevectors/. ./examples/toys_dataset/symb_statevectors - name: Generate HTML docs env: FIREBASE_CERTIFICATE: ${{ secrets.FIREBASE_CERTIFICATE }} @@ -33,7 +39,6 @@ jobs: python -m pip install --upgrade pip apt-get -y install --fix-missing git-core apt-get -y install build-essential - pip install -r doc/requirements.txt - name: Upload generated HTML as artifact uses: actions/upload-artifact@v4 with: diff --git a/README.md b/README.md index 785e2728..537436a0 100644 --- a/README.md +++ b/README.md @@ -123,12 +123,21 @@ To check the installation, open a python shell and type: import pyriemann_qiskit ``` -To enable Qiskit GPU optimization when using quantum simulation, run: +To enable Qiskit GPU optimization (for Linux) when using quantum simulation, run: + +``` +pip install .[optim_linux] +``` + +To use symbolic quantum simulation, run: ``` pip install .[optim] ``` +Which will enable [qiskit-symb](https://github.com/SimoneGasperini/qiskit-symb) +integration. + Note, Qiskit only provide binaries for Linux. For other platforms, or if you want to enable specific NVIDIA optimization for quantum, you need to build the binary [yourself](https://github.com/Qiskit/qiskit-aer/blob/main/CONTRIBUTING.md#building-with-gpu-support). diff --git a/doc/requirements.txt b/doc/requirements.txt index 26d6c462..cea7bf98 100644 --- a/doc/requirements.txt +++ b/doc/requirements.txt @@ -23,3 +23,5 @@ moabb==1.1.0 pyriemann==0.6 docplex>=2.21.207 firebase_admin==6.5.0 +qiskit-symb +symengine==0.11.0 diff --git a/examples/ERP/plot_classify_P300_bi.py b/examples/ERP/plot_classify_P300_bi.py index a8c70930..78616e49 100644 --- a/examples/ERP/plot_classify_P300_bi.py +++ b/examples/ERP/plot_classify_P300_bi.py @@ -86,7 +86,10 @@ # the non-qunatum SVM version used in qiskit # On a real Quantum computer (n_components = qubits) dim_red=PCA(n_components=5), - # params={'q_account_token': ''} + params={ + "n_jobs": 1, # Number of jobs for the simulator + # 'q_account_token': '' + }, ) # Here we provide a pipeline for comparison: diff --git a/examples/other_datasets/plot_financial_data.py b/examples/other_datasets/plot_financial_data.py index 9c75f36c..54a01a04 100644 --- a/examples/other_datasets/plot_financial_data.py +++ b/examples/other_datasets/plot_financial_data.py @@ -50,6 +50,7 @@ from sklearn.metrics import balanced_accuracy_score from pyriemann_qiskit.classification import QuanticSVM +from pyriemann_qiskit.utils.hyper_params_factory import gen_zz_feature_map from pyriemann_qiskit.utils.preprocessing import NdRobustScaler print(__doc__) @@ -370,7 +371,15 @@ def transform(self, X): # for the quantum SVM as for the classical one. gs.best_estimator_.steps[-1] = ( "quanticsvm", - QuanticSVM(quantum=True, C=best_C, gamma=best_gamma, seed=42), + QuanticSVM( + quantum=True, + C=best_C, + gamma=best_gamma, + gen_feature_map=gen_zz_feature_map(), + seed=42, + n_jobs=1, + use_qiskit_symb=False, + ), ) train_pred_qsvm = gs.best_estimator_.fit(X_train, y_train).predict(X_train) train_score_qsvm = balanced_accuracy_score(y_train, train_pred_qsvm) diff --git a/pyriemann_qiskit/classification.py b/pyriemann_qiskit/classification.py index 42da0e10..9caf9502 100644 --- a/pyriemann_qiskit/classification.py +++ b/pyriemann_qiskit/classification.py @@ -83,7 +83,7 @@ class QuanticClassifierBase(BaseEstimator, ClassifierMixin): If true, will output all intermediate results and logs. shots : int, default=1024 Number of repetitions of each circuit, for sampling. - gen_feature_map : Callable[int, QuantumCircuit | FeatureMap], \ + gen_feature_map : Callable[[int, str], QuantumCircuit | FeatureMap], \ default=Callable[int, ZZFeatureMap] Function generating a feature map to encode data into a quantum state. seed : int | None, default=None @@ -328,6 +328,9 @@ class QuanticSVM(QuanticClassifierBase): Predict is now using predict_proba with a softmax, when using QSVC. .. versionchanged:: 0.3.0 Add use_fidelity_state_vector_kernel parameter + .. versionchanged:: 0.4.0 + Add n_jobs and use_qiskit_symb parameter + for SymbFidelityStatevectorKernel Parameters ---------- @@ -359,17 +362,24 @@ class QuanticSVM(QuanticClassifierBase): If true, will output all intermediate results and logs. shots : int, default=1024 Number of repetitions of each circuit, for sampling. - gen_feature_map : Callable[int, QuantumCircuit | FeatureMap], \ + gen_feature_map : Callable[[int, str], QuantumCircuit | FeatureMap], \ default=Callable[int, ZZFeatureMap] Function generating a feature map to encode data into a quantum state. seed : int | None, default=None Random seed for the simulation - use_fidelity_state_vector_kernel: boolean (default=True) + use_fidelity_state_vector_kernel: boolean, default=True if True, use a FidelitystatevectorKernel for simulation. + use_qiskit_symb: boolean, default=True + This flag is used only if qiskit-symb is installed, and pegasos is False. + If True and the number of qubits < 9, then qiskit_symb is used. + n_jobs: boolean + The number of jobs for the qiskit-symb fidelity state vector + (if applicable) See Also -------- QuanticClassifierBase + SymbFidelityStatevectorKernel References ---------- @@ -407,6 +417,8 @@ def __init__( gen_feature_map=gen_zz_feature_map(), seed=None, use_fidelity_state_vector_kernel=True, + use_qiskit_symb=True, + n_jobs=4, ): QuanticClassifierBase.__init__( self, quantum, q_account_token, verbose, shots, gen_feature_map, seed @@ -416,14 +428,19 @@ def __init__( self.max_iter = max_iter self.pegasos = pegasos self.use_fidelity_state_vector_kernel = use_fidelity_state_vector_kernel + self.n_jobs = n_jobs + self.use_qiskit_symb = use_qiskit_symb def _init_algo(self, n_features): self._log("SVM initiating algorithm") if self.quantum: quantum_kernel = get_quantum_kernel( self._feature_map, + self.gen_feature_map, self._quantum_instance, self.use_fidelity_state_vector_kernel, + self.use_qiskit_symb and not self.pegasos, + self.n_jobs, ) if self.pegasos: self._log("[Warning] `gamma` is not supported by PegasosQSVC") @@ -498,7 +515,7 @@ class QuanticVQC(QuanticClassifierBase): If true, will output all intermediate results and logs shots : int, default=1024 Number of repetitions of each circuit, for sampling - gen_feature_map : Callable[int, QuantumCircuit | FeatureMap], \ + gen_feature_map : Callable[[int, str], QuantumCircuit | FeatureMap], \ default=Callable[int, ZZFeatureMap] Function generating a feature map to encode data into a quantum state. seed : int | None, default=None diff --git a/pyriemann_qiskit/pipelines.py b/pyriemann_qiskit/pipelines.py index d643fa8e..ea994300 100644 --- a/pyriemann_qiskit/pipelines.py +++ b/pyriemann_qiskit/pipelines.py @@ -458,7 +458,7 @@ class QuantumMDMVotingClassifier(BasePipeline): If true, will output all intermediate results and logs. shots : int (default:1024) Number of repetitions of each circuit, for sampling. - gen_feature_map : Callable[int, QuantumCircuit | FeatureMap] \ + gen_feature_map : Callable[[int, str], QuantumCircuit | FeatureMap] \ (default : Callable[int, ZZFeatureMap]) Function generating a feature map to encode data into a quantum state. diff --git a/pyriemann_qiskit/utils/hyper_params_factory.py b/pyriemann_qiskit/utils/hyper_params_factory.py index ec5140ae..c0f78052 100644 --- a/pyriemann_qiskit/utils/hyper_params_factory.py +++ b/pyriemann_qiskit/utils/hyper_params_factory.py @@ -24,8 +24,12 @@ def gen_x_feature_map(reps=2): Returns ------- - ret : XFeatureMap - An instance of XFeatureMap. + ret : Callable[[int, str], XFeatureMap] + A callable that takes into arguments: + + - the number of features + - the prefix string for the parameters + And returns an instance of XFeatureMap. Raises ------ @@ -40,16 +44,18 @@ def gen_x_feature_map(reps=2): Notes ----- .. versionadded:: 0.2.0 + .. versionchanged:: 0.4.0 + Possibility to specify parameter prefix. """ if reps < 1: raise ValueError(f"Parameter reps must be superior or equal to 1 (Got {reps})") - return lambda n_features: PauliFeatureMap( + return lambda n_features, param_prefix="x": PauliFeatureMap( feature_dimension=n_features, paulis=["X"], reps=reps, data_map_func=None, - parameter_prefix="x", + parameter_prefix=param_prefix, insert_barriers=False, name="XFeatureMap", ) @@ -69,8 +75,12 @@ def gen_z_feature_map(reps=2): Returns ------- - ret : ZFeatureMap - An instance of ZFeatureMap. + ret : Callable[[int, str], ZFeatureMap] + A callable that takes into arguments: + + - the number of features + - the prefix string for the parameters + And returns an instance of ZFeatureMap. Raises ------ @@ -85,11 +95,17 @@ def gen_z_feature_map(reps=2): Notes ----- .. versionadded:: 0.2.0 + .. versionchanged:: 0.4.0 + Possibility to specify parameter prefix. """ if reps < 1: raise ValueError(f"Parameter reps must be superior or equal to 1 (Got {reps})") - return lambda n_features: ZFeatureMap(feature_dimension=n_features, reps=reps) + return lambda n_features, param_prefix="x": ZFeatureMap( + feature_dimension=n_features, + reps=reps, + parameter_prefix=param_prefix, + ) def gen_zz_feature_map(reps=2, entanglement="linear"): @@ -113,8 +129,12 @@ def gen_zz_feature_map(reps=2, entanglement="linear"): Returns ------- - ret : ZZFeatureMap - An instance of ZZFeatureMap. + ret : Callable[[int, str], ZZFeatureMap] + A callable that takes into arguments: + + - the number of features + - the prefix string for the parameters + And returns an instance of ZZFeatureMap. Raises ------ @@ -131,12 +151,17 @@ def gen_zz_feature_map(reps=2, entanglement="linear"): Notes ----- .. versionadded:: 0.0.1 + .. versionchanged:: 0.4.0 + Possibility to specify parameter prefix. """ if reps < 1: raise ValueError(f"Parameter reps must be superior or equal to 1 (Got {reps})") - return lambda n_features: ZZFeatureMap( - feature_dimension=n_features, reps=reps, entanglement=entanglement + return lambda n_features, param_prefix="x": ZZFeatureMap( + feature_dimension=n_features, + reps=reps, + entanglement=entanglement, + parameter_prefix=param_prefix, ) diff --git a/pyriemann_qiskit/utils/quantum_provider.py b/pyriemann_qiskit/utils/quantum_provider.py index 10e1fe9f..de7f787a 100644 --- a/pyriemann_qiskit/utils/quantum_provider.py +++ b/pyriemann_qiskit/utils/quantum_provider.py @@ -1,8 +1,12 @@ """Module containing helpers for IBM quantum backends providers and simulators.""" +import joblib import logging +import os +import pickle +import numpy as np from qiskit_aer import AerSimulator from qiskit_aer.quantum_info import AerStatevector from qiskit_algorithms.state_fidelities import ComputeUncompute @@ -12,6 +16,122 @@ FidelityQuantumKernel, ) +try: + from qiskit_symb.quantum_info import Statevector + + QISKIT_SYMB = True +except ImportError: + QISKIT_SYMB = False + + +class SymbFidelityStatevectorKernel: + + """Symbolic Statevector kernel + + An implementation of the quantum kernel for classically simulated + state vectors [1]_ using qiskit-symb for symbolic representation + of statevectors [2]_. + + Here, the kernel function is defined as the overlap of two simulated quantum + statevectors produced by a parametrized quantum circuit (called feature map) [1]_. + + Notes + ----- + .. versionadded:: 0.4.0 + + Parameters + ---------- + feature_map: QuantumCircuit | FeatureMap + An instance of a feature map. + gen_feature_map: Callable[[int, str], QuantumCircuit | FeatureMap], \ + default=Callable[int, ZZFeatureMap] + Function generating a feature map to encode data into a quantum state. + n_jobs: int + The number of job for fidelity evaluation. + If null or negative, the number of jobs is set to 1 + If set to 1, evaluation will run on the main thread. + + References + ---------- + .. [1] \ + https://github.com/qiskit-community/qiskit-machine-learning/blob/30dad803e9457f955464220eddc1e55a65452bbc/qiskit_machine_learning/kernels/fidelity_statevector_kernel.py#L31 + .. [2] https://github.com/SimoneGasperini/qiskit-symb/issues/6 + + """ + + def __init__(self, feature_map, gen_feature_map, n_jobs=1): + self.n_jobs = n_jobs if n_jobs >= 1 else 1 + cached_file = os.path.join( + "symb_statevectors", f"{feature_map.name}-{feature_map.reps}" + ) + + if os.path.isfile(cached_file): + print("Loading symbolic Statevector from cache") + file = open(cached_file, "rb") + sv = pickle.load(file) + else: + print("Computing symbolic Statevector") + fm2 = gen_feature_map(feature_map.num_qubits, "b") + self.circuit = feature_map.compose(fm2.inverse()).decompose() + sv = Statevector(self.circuit) + print(f"Dumping to {cached_file}") + file = open(cached_file, "wb") + pickle.dump(sv, file) + + self.function = sv.to_lambda() + + def evaluate(self, x_vec, y_vec=None): + """Evaluate the quantum kernel. + + Returns + ------- + kernel : ndarray, shape (len(x_vec), len(y_vec)) + The kernel matrix. + + Notes + ----- + .. versionadded:: 0.4.0 + """ + if y_vec is None: + y_vec = x_vec + + x_vec_len = len(x_vec) + y_vec_len = len(y_vec) + + is_sim = x_vec_len == y_vec_len and (x_vec == y_vec).all() + + kernel_matrix = np.zeros((x_vec_len, y_vec_len)) + + chunck = x_vec_len // self.n_jobs + + def compute_fidelity_partial_matrix(i_thread): + for i in range(i_thread * chunck, (i_thread + 1) * chunck): + x = x_vec[i] + for j in range(i if is_sim else y_vec_len): + y = y_vec[j] + if isinstance(x, np.float64): + # Pegagos implementation + fidelity = abs(self.function(x, y)[0, 0]) ** 2 + else: + fidelity = abs(self.function(*x, *y)[0, 0]) ** 2 + + kernel_matrix[i, j] = fidelity + if is_sim: + kernel_matrix[j, i] = fidelity + return kernel_matrix + + if self.n_jobs == 1: + return compute_fidelity_partial_matrix(0) + else: + print("n_jobs greater than 1, parallelizing") + results = joblib.Parallel(n_jobs=self.n_jobs)( + joblib.delayed(compute_fidelity_partial_matrix)(i_thread) + for i_thread in range(self.n_jobs) + ) + for result in results: + kernel_matrix += result + return kernel_matrix + def get_provider(): """Return an IBM quantum provider. @@ -93,12 +213,22 @@ def get_device(provider, min_qubits): ) -def get_quantum_kernel(feature_map, quantum_instance, use_fidelity_state_vector_kernel): +def get_quantum_kernel( + feature_map, + gen_feature_map, + quantum_instance, + use_fidelity_state_vector_kernel, + use_qiskit_symb, + n_jobs=4, +): """Get a quantum kernel Return an instance of FidelityQuantumKernel or FidelityStatevectorKernel (in the case of a simulation). + For simulation with a small number of qubits (< 9), and `use_qiskit_symb` is True, + qiskit-symb is used. + Parameters ---------- feature_map: QuantumCircuit | FeatureMap @@ -106,36 +236,60 @@ def get_quantum_kernel(feature_map, quantum_instance, use_fidelity_state_vector_ quantum_instance: BaseSampler A instance of BaseSampler. use_fidelity_state_vector_kernel: boolean - if True, use a FidelitystatevectorKernel for simulation. + If True, use a FidelitystatevectorKernel for simulation. + use_qiskit_symb: boolean + This flag is used only if qiskit-symb is installed. + If True and the number of qubits < 9, then qiskit_symb is used. + n_jobs: boolean + The number of jobs for the qiskit-symb fidelity state vector + (if applicable) Returns ------- kernel: QuantumKernel The quantum kernel. + See also + -------- + SymbFidelityStatevectorKernel + Notes ----- .. versionadded:: 0.3.0 + .. versionchanged:: 0.4.0 + Add support for qiskit-symb """ if use_fidelity_state_vector_kernel and isinstance( quantum_instance._backend, AerSimulator ): - logging.log( - logging.WARN, - """FidelityQuantumKernel skipped because of time. + # For simulation: + if QISKIT_SYMB and feature_map.num_qubits <= 9 and use_qiskit_symb: + # With a small number of qubits, let's use qiskit-symb + # See: + # https://medium.com/qiskit/qiskit-symb-a-qiskit-ecosystem-package-for-symbolic-quantum-computation-b6b4407fa705 + kernel = SymbFidelityStatevectorKernel( + feature_map, gen_feature_map, n_jobs=n_jobs + ) + logging.log( + logging.WARN, + """Using SymbFidelityStatevectorKernel""", + ) + else: + # For a larger number of qubits, + # we will not use FidelityQuantumKernel as it is slow. See + # https://github.com/qiskit-community/qiskit-machine-learning/issues/547#issuecomment-1486527297 + kernel = FidelityStatevectorKernel( + feature_map=feature_map, + statevector_type=AerStatevector, + shots=quantum_instance.options["shots"], + ) + logging.log( + logging.WARN, + """FidelityQuantumKernel skipped because of time. Using FidelityStatevectorKernel with AerStatevector. Seed cannot be set with FidelityStatevectorKernel. Increase the number of shots to diminish the noise.""", - ) - - # if this is a simulation, - # we will not use FidelityQuantumKernel as it is slow. See - # https://github.com/qiskit-community/qiskit-machine-learning/issues/547#issuecomment-1486527297 - kernel = FidelityStatevectorKernel( - feature_map=feature_map, - statevector_type=AerStatevector, - shots=quantum_instance.options["shots"], - ) + ) else: kernel = FidelityQuantumKernel( feature_map=feature_map, fidelity=ComputeUncompute(quantum_instance) diff --git a/setup.py b/setup.py index e636d7f0..243bce8c 100644 --- a/setup.py +++ b/setup.py @@ -49,7 +49,7 @@ 'firebase_admin==6.5.0', 'scikit-learn==1.5.1', 'tqdm', - 'pandas' + 'pandas', ], extras_require={'docs': [ 'sphinx-gallery', @@ -64,6 +64,7 @@ 'tests': ['pytest', 'seaborn', 'flake8', 'mne', 'pooch'], # GPU optimization not available on all platform. # See https://github.com/Qiskit/qiskit-aer/issues/929#issuecomment-691716936 - 'optim': ['qiskit-aer-gpu==0.15.0']}, + 'optim': ['qiskit-symb', 'symengine==0.11.0'], + 'optim_linux': ['qiskit-aer-gpu==0.15.0']}, zip_safe=False, ) diff --git a/symb_statevectors/XFeatureMap-2 b/symb_statevectors/XFeatureMap-2 new file mode 100644 index 00000000..d4a034e3 Binary files /dev/null and b/symb_statevectors/XFeatureMap-2 differ diff --git a/symb_statevectors/XFeatureMap-3 b/symb_statevectors/XFeatureMap-3 new file mode 100644 index 00000000..1f662c74 Binary files /dev/null and b/symb_statevectors/XFeatureMap-3 differ diff --git a/symb_statevectors/ZFeatureMap-2 b/symb_statevectors/ZFeatureMap-2 new file mode 100644 index 00000000..2c388627 Binary files /dev/null and b/symb_statevectors/ZFeatureMap-2 differ diff --git a/symb_statevectors/ZFeatureMap-3 b/symb_statevectors/ZFeatureMap-3 new file mode 100644 index 00000000..90296404 Binary files /dev/null and b/symb_statevectors/ZFeatureMap-3 differ diff --git a/symb_statevectors/ZZFeatureMap-2 b/symb_statevectors/ZZFeatureMap-2 new file mode 100644 index 00000000..085b8943 Binary files /dev/null and b/symb_statevectors/ZZFeatureMap-2 differ diff --git a/symb_statevectors/ZZFeatureMap-3 b/symb_statevectors/ZZFeatureMap-3 new file mode 100644 index 00000000..c0590a04 Binary files /dev/null and b/symb_statevectors/ZZFeatureMap-3 differ