Repository: quantumlib/OpenFermion
Pull request #558 (2020-02-26 20:19:26)

Prony estimator
A slimmed down version of Prony's method for QPE signal processing.

obriente: Seems that this is a case of the bug in the coverage script.

Patch:
diff --git a/src/openfermion/measurements/__init__.py b/src/openfermion/measurements/__init__.py
index 84a4d58..60c67a3 100644
--- a/src/openfermion/measurements/__init__.py
+++ b/src/openfermion/measurements/__init__.py
@@ -21,3 +21,5 @@ from ._rdm_equality_constraints import (one_body_fermion_constraints,
 from ._qubit_partitioning import (binary_partition_iterator,
                                   partition_iterator,
                                   pauli_string_iterator)
+
+from ._prony import (prony)
diff --git a/src/openfermion/measurements/_prony.py b/src/openfermion/measurements/_prony.py
new file mode 100644
index 0000000..8066098
--- /dev/null
+++ b/src/openfermion/measurements/_prony.py
@@ -0,0 +1,55 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+"""An implementation of Prony's method (or the matrix pencil method)
+This fits a signal f(t) to sum_i=1^M a_i gamma_i^t, where a_i, gamma_i
+are complex numbers
+"""
+import scipy
+import numpy
+
+
+def prony(signal):
+    """Estimates amplitudes and phases of a sparse signal using Prony's method.
+
+    Single-ancilla quantum phase estimation returns a signal
+    g(k)=sum (aj*exp(i*k*phij)), where aj and phij are the amplitudes
+    and corresponding eigenvalues of the unitary whose phases we wish
+    to estimate. When more than one amplitude is involved, Prony's method
+    provides a simple estimation tool, which achieves near-Heisenberg-limited
+    scaling (error scaling as N^{-1/2}K^{-3/2}).
+
+    Args:
+        signal(1d complex array): the signal to fit
+
+    Returns:
+        amplitudes(list of complex values): the amplitudes a_i, in descending order
+            by their complex magnitude
+        phases(list of complex values): the complex frequencies gamma_i,
+            correlated with amplitudes.
+    """
+
+    num_freqs = len(signal) // 2
+    hankel0 = scipy.linalg.hankel(c=signal[:num_freqs],
+                                  r=signal[num_freqs - 1:-1])
+    hankel1 = scipy.linalg.hankel(c=signal[1:num_freqs + 1],
+                                  r=signal[num_freqs:])
+    shift_matrix = scipy.linalg.lstsq(hankel0.T, hankel1.T)[0]
+    phases = numpy.linalg.eigvals(shift_matrix.T)
+
+    generation_matrix = numpy.array(
+        [[phase**k for phase in phases] for k in range(len(signal))])
+    amplitudes = scipy.linalg.lstsq(generation_matrix, signal)[0]
+
+    amplitudes, phases = zip(*sorted(
+        zip(amplitudes, phases), key=lambda x: numpy.abs(x[0]), reverse=True))
+
+    return numpy.array(amplitudes), numpy.array(phases)


Test patch:
diff --git a/src/openfermion/measurements/_prony_test.py b/src/openfermion/measurements/_prony_test.py
new file mode 100644
index 0000000..9b2dacc
--- /dev/null
+++ b/src/openfermion/measurements/_prony_test.py
@@ -0,0 +1,49 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+""" Tests for _prony.py"""
+import pytest
+from ._prony import prony
+import numpy
+
+
+def test_zeros():
+    signal = numpy.zeros(10)
+    amplitudes, phases = prony(signal)
+    assert (len(amplitudes) == 5)
+    assert (len(phases) == 5)
+    for j in range(5):
+        numpy.testing.assert_allclose(amplitudes[j], 0)
+        numpy.testing.assert_allclose(phases[j], 0)
+
+
+def test_signal():
+    x_vec = numpy.linspace(0, 1, 11)
+    y_vec = (0.5 * numpy.exp(1j * x_vec * 3) + 0.3 * numpy.exp(1j * x_vec * 5) +
+             0.15 * numpy.exp(1j * x_vec * 1.5) +
+             0.1 * numpy.exp(1j * x_vec * 4) +
+             0.05 * numpy.exp(1j * x_vec * 1.2))
+    print(y_vec)
+    amplitudes, phases = prony(y_vec)
+    assert (len(amplitudes) == 5)
+    assert (len(phases) == 5)
+    for a, p in zip(amplitudes, phases):
+        print(a, numpy.angle(p))
+    numpy.testing.assert_allclose(numpy.abs(amplitudes[0]), 0.5, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.abs(amplitudes[1]), 0.3, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.abs(amplitudes[2]), 0.15, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.abs(amplitudes[3]), 0.1, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.abs(amplitudes[4]), 0.05, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.angle(phases[0]), 0.3, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.angle(phases[1]), 0.5, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.angle(phases[2]), 0.15, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.angle(phases[3]), 0.4, atol=1e-6)
+    numpy.testing.assert_allclose(numpy.angle(phases[4]), 0.12, atol=1e-6)