【代码贡献】Fix QPCA crash on 4-feature input and k=2 returning 0

pyqpanda-algorithm/pyqpanda_alg/QPCA/QPCA.py

@@ -13,7 +13,6 @@
 from pyqpanda3.core import CPUQVM, QCircuit, QProg, TOFFOLI, SWAP, CNOT, U1, U3, CR, I, H, X, Y, measure
 import numpy as np
 import math
-import sys
 
 from .. plugin import *
 
@@ -201,15 +200,22 @@ def qpca(sample_A, k):
     """
     QPCA is a quantum version of the classical PCA algorithm, which is widely used in data analysis and machine learning.
 
+    Only 2-feature input is supported. The n=4 branches of the circuit helpers
+    are incomplete, so a larger input raises NotImplementedError.
+
     Parameters:
         sample_A: ``ndarray``\n
-            the input matrix for analysis
+            the input matrix for analysis, shape (n_samples, 2)
         k: ``int``\n
-            the dimension to reduce
+            the dimension to reduce, 1 or 2
 
     Returns:
         out: ``ndarray``\n
-            the output matrix after reducing dimension
+            the output matrix after reducing dimension, shape (n_samples, k)
+
+    Raises:
+        ValueError: k is not 1 or 2.\n
+        NotImplementedError: sample_A has other than 2 features.
 
     Examples:
         .. code-block:: python
@@ -222,61 +228,55 @@ def qpca(sample_A, k):
             print(data_q)
     """
     norm_x, A = _preprocessing(sample_A)
-    lambda_A, vector_A = np.linalg.eig(A)
-    A1 = A.reshape(1, 2 ** A.shape[0])[0]
+    n = A.shape[0]
+    if n != 2:
+        raise NotImplementedError(
+            'qpca supports 2-feature input only, got %d features' % n)
+    if k not in (1, 2):
+        raise ValueError('k must be 1 or 2, got %r' % (k,))
+
+    # the covariance matrix is symmetric, so eigh gives real eigenvalues and an
+    # orthonormal basis. sorted descending, principal component first.
+    lambda_A, vector_A = np.linalg.eigh(A)
+    order = np.argsort(lambda_A)[::-1]
+    lambda_A = lambda_A[order]
+    vector_A = vector_A[:, order]
+    A1 = A.reshape(1, 2 ** n)[0]
     state_vector = np.sqrt(1 / np.sum(A1 * A1)) * A1
 
-    if A.shape[0] == 2:
-        if k == 2:
-            tao = min(lambda_A) - 1
-        elif k == 1:
-            tao = np.sum(lambda_A)/2
-        else:
-            print('The K Error!')
-            sys.exit()
-        qm = _QMachine(5, 5)
-    if A.shape[0] == 4:
-        qm = _QMachine(8, 8)
+    # tao is the eigenvalue threshold. k=1 keeps the top eigenvalue only,
+    # k=2 sits below both and keeps the whole space.
+    if k == 2:
+        tao = min(lambda_A) - 1
+    else:
+        tao = np.sum(lambda_A) / 2
+    qm = _QMachine(5, 5)
 
     prog = QProg()
     cir = QCircuit()
-    _init_cir(qm, state_vector, A.shape[0])
-    cir << _phase_estimation_cir(qm.q_list, lambda_A, tao, A.shape[0])
-    cir << _transition_cir(qm.q_list, A.shape[0])
+    _init_cir(qm, state_vector, n)
+    cir << _phase_estimation_cir(qm.q_list, lambda_A, tao, n)
+    cir << _transition_cir(qm.q_list, n)
     cir << _cnot_cir(qm.q_list)
-    cir << _transition_reverse_cir(qm.q_list, A.shape[0])
-    cir << _phase_estimation_reverse_cir(qm.q_list, lambda_A, tao, A.shape[0])
+    cir << _transition_reverse_cir(qm.q_list, n)
+    cir << _phase_estimation_reverse_cir(qm.q_list, lambda_A, tao, n)
     prog << cir
     prog <<_measure_cir(prog, qm.q_list, qm.q_list)
-    # result = qm.machine.run_with_configuration(prog, qm.c_list, 8192)
     qm.machine.run(prog, 8192)
     result = qm.machine.result().get_prob_dict(qm.q_list)
-    a = []
-    data = 0
-    if A.shape[0] == 2:
-        for i, v in enumerate(result.keys()):
-            if int(v[-1]) == 1:
-                # a.append(float("%.4f" % (result[v] / 8192)))
-                a.append(float("%.4f" % (result[v])))
-        if k == 2:
-            result_idealA = state_vector
-        if k == 1:
-            i = np.argmax(lambda_A)
-            result_idealA = (np.kron(vector_A[:, i], vector_A[:, i]) * lambda_A[i]) / (
-                np.sqrt(lambda_A[i] * lambda_A[i]))
-        result_idealA1 = []
-        for i in range(len(result_idealA)):
-            if result_idealA[i] != 0:
-                result_idealA1.append(result_idealA[i])
-        result_idealA = result_idealA1
-        sum_A = np.sum(np.array(a))
-        result_circuitA = []
-        for i in range(len(a)):
-            result_circuitA.append(np.sqrt(a[i] / sum_A))
-        if k == 1:
-            vector_qpca = []
-            vector_qpca.append(np.sqrt(result_circuitA[0]))
-            vector_qpca.append(np.sqrt(result_circuitA[-1]))
-            vector_qpca = np.array(vector_qpca).reshape(1, 2)
-            data = np.dot(norm_x, np.transpose(vector_qpca))
-    return data
+
+    if k == 2:
+        # tao is below every eigenvalue, so post-selection succeeds with
+        # probability 1 and nothing is filtered out. the retained subspace is
+        # the whole feature space and its basis comes from vector_A.
+        return np.dot(norm_x, vector_A)
+
+    # post-select the ancilla, qubit 0, which is the last character of the key.
+    # sorted() so the branch order does not depend on dict insertion order.
+    a = [result[v] for v in sorted(result) if int(v[-1]) == 1]
+    sum_A = np.sum(np.array(a))
+    # probabilities back to amplitudes. unit norm when post-selection leaves
+    # two outcomes, which is the non-degenerate case.
+    result_circuitA = [np.sqrt(x / sum_A) for x in a]
+    vector_qpca = np.array([result_circuitA[0], result_circuitA[-1]]).reshape(1, 2)
+    return np.dot(norm_x, np.transpose(vector_qpca))

test/QPCA/Test_qpca.py

@@ -8,6 +8,15 @@
 sys.path.append((Path.cwd().parent.parent).__str__())
 
 
+def _classical_pca(x):
+    """Reference PCA: eigendecomposition of the covariance matrix, descending."""
+    norm_x = x - x.mean(axis=0)
+    cov = np.dot(norm_x.T, norm_x)
+    w, v = np.linalg.eigh(cov)
+    order = np.argsort(w)[::-1]
+    return norm_x, w[order], v[:, order]
+
+
 class TestQPCA:
     """QPCA测试类"""
 
@@ -16,6 +25,11 @@ def sample_data(self):
         """提供标准测试数据"""
         return np.array([[-1, 2], [-2, -1], [-1, -2], [1, 3], [2, 1], [3, 2]])
 
+    @pytest.fixture
+    def sample_data_4d(self):
+        rng = np.random.RandomState(0)
+        return rng.randn(20, 4)
+
     def test_qpca_normal(self, sample_data):
         """测试QPCA正常功能"""
         from pyqpanda_alg.QPCA import qpca
@@ -26,5 +40,64 @@ def test_qpca_normal(self, sample_data):
         assert data_q.shape == (6,1)
         assert isinstance(data_q, np.ndarray)
 
+    def test_qpca_k1_matches_classical_pca(self, sample_data):
+        from pyqpanda_alg.QPCA import qpca
+
+        data_q = qpca(sample_data, 1)[:, 0]
+        norm_x, _, v = _classical_pca(sample_data)
+        expected = np.dot(norm_x, v[:, 0])
+        if np.dot(data_q, expected) < 0:
+            expected = -expected
+
+        # sampling noise on 8192 shots keeps this loose
+        assert np.corrcoef(data_q, expected)[0, 1] > 0.999
+        assert np.max(np.abs(data_q - expected)) < 0.15
+
+    def test_qpca_k1_direction_is_unit_norm(self, sample_data):
+        """Holds for this fixture, not in general: a degenerate spectrum leaves
+        the post-selection with a single outcome and the norm comes out at
+        sqrt(2). Pinned to the documented example on purpose."""
+        from pyqpanda_alg.QPCA import qpca
+
+        norm_x, _, _ = _classical_pca(sample_data)
+        data_q = qpca(sample_data, 1)
+        # recover the projection direction the circuit produced
+        vec = np.linalg.lstsq(norm_x, data_q, rcond=None)[0][:, 0]
+        assert np.linalg.norm(vec) == pytest.approx(1.0, abs=1e-6)
+
+    def test_qpca_k2_returns_pca_scores(self, sample_data):
+        from pyqpanda_alg.QPCA import qpca
+
+        data_q = qpca(sample_data, 2)
+        norm_x, w, v = _classical_pca(sample_data)
+
+        assert isinstance(data_q, np.ndarray)
+        assert data_q.shape == (6, 2)
+
+        # PCA scores are decorrelated and carry the eigenvalues as their
+        # column variances, in descending order. neither holds for the
+        # centred input, so this pins the basis and not just the subspace.
+        cov_q = np.dot(data_q.T, data_q)
+        assert abs(cov_q[0, 1]) < 1e-9
+        assert np.allclose(np.diag(cov_q), w)
+
+        expected = np.dot(norm_x, v)
+        for j in range(2):
+            sign = np.sign(np.dot(data_q[:, j], expected[:, j]))
+            assert np.allclose(data_q[:, j], sign * expected[:, j])
+
+    @pytest.mark.parametrize('k', [0, 3, -1])
+    def test_qpca_invalid_k_raises(self, sample_data, k):
+        from pyqpanda_alg.QPCA import qpca
+
+        with pytest.raises(ValueError):
+            qpca(sample_data, k)
+
+    def test_qpca_four_features_raises(self, sample_data_4d):
+        from pyqpanda_alg.QPCA import qpca
 
+        with pytest.raises(NotImplementedError):
+            qpca(sample_data_4d, 1)
 
+        with pytest.raises(NotImplementedError):
+            qpca(sample_data_4d, 2)