Fix QSVM problem wit non-positive semi-definite kernel matrices (qisk…

…it-community/qiskit-aqua#1190)

* Update _qsvm_binary.py

* move psd approx to static method for kernel matrix creation

* only apply PSD approximation for symmetric matrices

* remove `print(mat)` from qsvm.py `get_kernel_matrix`

* Add unit test to check enforce_psd

Co-authored-by: Manoel Marques <Manoel.Marques@ibm.com>
Co-authored-by: woodsp <woodsp@us.ibm.com>
Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>
Co-authored-by: Julien Gacon <gaconju@gmail.com>

qiskit/aqua/algorithms/classifiers/qsvm/qsvm.py

@@ -215,7 +215,7 @@ def construct_circuit(self, x1, x2, measurement=False):
         return QSVM._construct_circuit((x1, x2), self.feature_map, measurement)
 
     @staticmethod
-    def get_kernel_matrix(quantum_instance, feature_map, x1_vec, x2_vec=None):
+    def get_kernel_matrix(quantum_instance, feature_map, x1_vec, x2_vec=None, enforce_psd=True):
         """
         Construct kernel matrix, if x2_vec is None, self-innerproduct is conducted.
 
@@ -234,6 +234,9 @@ def get_kernel_matrix(quantum_instance, feature_map, x1_vec, x2_vec=None):
                                     D is the feature dimension
             x2_vec (numpy.ndarray): data points, 2-D array, N2xD, where N2 is the number of data,
                                     D is the feature dimension
+            enforce_psd (bool): enforces that the kernel matrix is positive semi-definite by setting
+                                negative eigenvalues to zero. This is only applied in the symmetric
+                                case, i.e., if `x2_vec == None`.
         Returns:
             numpy.ndarray: 2-D matrix, N1xN2
         """
@@ -357,6 +360,14 @@ def get_kernel_matrix(quantum_instance, feature_map, x1_vec, x2_vec=None):
                     if is_symmetric:
                         mat[j, i] = mat[i, j]
 
+        if enforce_psd and is_symmetric and not is_statevector_sim:
+            # Find the closest positive semi-definite approximation to kernel matrix, in case it is
+            # symmetric. The (symmetric) matrix should always be positive semi-definite by
+            # construction, but this can be violated in case of noise, such as sampling noise, thus,
+            # the adjustment is only done if NOT using the statevector simulation.
+            D, U = np.linalg.eig(mat)
+            mat = U @ np.diag(np.maximum(0, D)) @ U.transpose()
+
         return mat
 
     def construct_kernel_matrix(self, x1_vec, x2_vec=None, quantum_instance=None):

test/aqua/test_qsvm.py

@@ -16,15 +16,19 @@
 
 import os
 from test.aqua import QiskitAquaTestCase
+
 import numpy as np
 from ddt import ddt, data
+
 from qiskit import BasicAer, QuantumCircuit
 from qiskit.circuit.library import ZZFeatureMap
 from qiskit.aqua import QuantumInstance, aqua_globals, MissingOptionalLibraryError
 from qiskit.aqua.components.multiclass_extensions import (ErrorCorrectingCode,
                                                           AllPairs,
                                                           OneAgainstRest)
 from qiskit.aqua.algorithms import QSVM
+from qiskit.aqua.utils import split_dataset_to_data_and_labels, optimize_svm
+from qiskit.ml.datasets import ad_hoc_data
 
 
 @ddt
@@ -218,3 +222,57 @@ def test_multiclass(self, multiclass_extension):
             self.assertEqual(result['predicted_classes'], predicted_classes[multiclass_extension])
         except MissingOptionalLibraryError as ex:
             self.skipTest(str(ex))
+
+    def test_matrix_psd(self):
+        """ Test kernel matrix positive semi-definite enforcement. """
+        try:
+            from cvxpy.error import DQCPError
+        except ImportError:
+            self.skipTest('cvxpy does not appeat to be installed')
+
+        seed = 10598
+        feature_dim = 2
+        _, training_input, _, _ = ad_hoc_data(
+            training_size=10,
+            test_size=5,
+            n=feature_dim,
+            gap=0.3
+        )
+        training_input, _ = split_dataset_to_data_and_labels(training_input)
+        training_data = training_input[0]
+        training_labels = training_input[1]
+        labels = training_labels * 2 - 1  # map label from 0 --> -1 and 1 --> 1
+        labels = labels.astype(np.float)
+
+        feature_map = ZZFeatureMap(feature_dimension=feature_dim, reps=2, entanglement='linear')
+
+        with self.assertRaises(DQCPError):
+            # Sampling noise means that the kernel matrix will not quite be positive
+            # semi-definite which will cause the optimize svm to fail
+            backend = BasicAer.get_backend('qasm_simulator')
+            quantum_instance = QuantumInstance(backend, shots=1024, seed_simulator=seed,
+                                               seed_transpiler=seed)
+            kernel_matrix = QSVM.get_kernel_matrix(quantum_instance, feature_map=feature_map,
+                                                   x1_vec=training_data, enforce_psd=False)
+            _ = optimize_svm(kernel_matrix, labels)
+
+        # This time we enforce that the matrix be positive semi-definite which runs logic to
+        # make it so.
+        backend = BasicAer.get_backend('qasm_simulator')
+        quantum_instance = QuantumInstance(backend, shots=1024, seed_simulator=seed,
+                                           seed_transpiler=seed)
+        kernel_matrix = QSVM.get_kernel_matrix(quantum_instance, feature_map=feature_map,
+                                               x1_vec=training_data, enforce_psd=True)
+        alpha, b, support = optimize_svm(kernel_matrix, labels)
+
+        expected_alpha = [0.855861781, 2.59807482, 0, 0.962959215,
+                          1.08141696, 0.217172547, 0, 0,
+                          0.786462904, 0, 0.969727949, 1.98066946,
+                          0, 0, 1.62049430, 0,
+                          0.394212728, 0, 0.507740935, 1.02910286]
+        expected_b = [-0.17543365]
+        expected_support = [True, True, False, True, True, True, False, False, True, False,
+                            True, True, False, False, True, False, True, False, True, True]
+        np.testing.assert_array_almost_equal(alpha, expected_alpha)
+        np.testing.assert_array_almost_equal(b, expected_b)
+        np.testing.assert_array_equal(support, expected_support)