This question follows the post https://quantumcomputing.stackexchange.com/posts/15070/.

I am working on an implementation of HHL algorithm, to do so, I need to map $|\lambda\rangle\mapsto |\arcsin(\frac{C}{\lambda})\rangle$, where $|\arcsin(\frac{C}{\lambda})\rangle$ is a binary representation $\arcsin(\frac{C}{\lambda})$ with $m$ qubits.

Recently, a circuit implementing a Piecewise Chebyshev approximation of a given function $f$ was implemented in Qiskit, documentation can be found here : https://qiskit.org/documentation/stubs/qiskit.circuit.library.PiecewiseChebyshev.html#qiskit-circuit-library-piecewisechebyshev.

I was working on the given example:

import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev

f_x, degree, breakpoints, num_state_qubits = lambda x: 0, 2, None, 2
pw_approximation = PiecewiseChebyshev(f_x, degree, breakpoints, num_state_qubits)
qc = QuantumCircuit(pw_approximation.num_qubits)
qc.append(pw_approximation.to_instruction(), qc.qubits)

which should implement a piecewise approximation of the zero function on two qubits ... We can take a look of what the circuit looks like:

enter image description here

the first and second qubits are state qubits and I guess that the three others are ancillary qubits.

The initial state of the system is set to $\frac{1}{2} (|00\rangle+|01\rangle+|10\rangle+|11\rangle)$ so we may look at all possible results of the approximation.

But while looking at measurements of the output qubits, I am facing weird results:

enter image description here

The third qubit is set to one, while the two first are unchanged. I tried a lot of different functions to understand how one can interpret the resulting qubits, but I still cannot figure out what the circuit does. Does anyone have an idea ?

  • $\begingroup$ Note: The bounty is assigned if you correctly explain (if and) how to use this component to implement the requested transformation |x⟩|0...0⟩→|x⟩|f(x)⟩ and NOT the transformation |x⟩|0⟩→|x⟩(cos(f(x))|0⟩+sin(f(x))|1⟩). The latter is being explained here. $\endgroup$
    – incud
    Jul 7 at 9:21
  • $\begingroup$ If you think that phase estimation is a necessary step to construct |f(x)>, please show how to put together PiecewiseChebyshev with PhaseEstimation which is already implemented in Qiskit $\endgroup$
    – incud
    Jul 7 at 9:23

The source of the problem was numpy.polynomial.chebyshev.interpolate, which is used to approximate the input function but raises problems when approximating a constant function.

The code in Qiskit has been changed adding a special case for such functions and explaining how constant functions should be declared, i.e. as f_x = constant rather than with the lambda declaration.

The updated code now performs as expected:

import numpy as np
from qiskit import *
from qiskit.visualization import plot_histogram
from qiskit.circuit.library.arithmetic.piecewise_chebyshev import PiecewiseChebyshev
f_x, degree, breakpoints, num_state_qubits = 0, 1, [0], 2
pw_approximation = PiecewiseChebyshev(f_x, degree, breakpoints, num_state_qubits)
qc = QuantumCircuit(pw_approximation.num_qubits, pw_approximation.num_qubits)
qc.append(pw_approximation.to_instruction(), qc.qubits)
qc.measure(list(range(pw_approximation.num_qubits)), list(range(pw_approximation.num_qubits)))

enter image description here

backend = BasicAer.get_backend('qasm_simulator') # the device to run on
result = execute(qc, backend, shots=10000).result()
counts  = result.get_counts(qc)

enter image description here

And here is the link to the PR in case you want to check its status.

Regarding the operation $f(x)|0>\mapsto |f(x)>$, it is not supported in Qiskit yet but here is the relevant literature.


Since the PiecewiseChebyshev object implements the transformation $U |x\rangle |0\rangle \to |x\rangle \Big(\cos(f(x))|0\rangle + \sin(f(x))|1\rangle\Big)$, you need to apply amplitude estimation process in order to retrieve the value $\sin(f(x))$ into some $m$-qubit quantum register. You can of course define $f$ such that the $\sin$ function disappear.

I have tried with $f = \arcsin \sqrt{x / 2^m}$ so that I can immediately check if the circuit works: if I encode as input the $n$-qubit state $|i\rangle$ then I will read in the output $m$-qubit quantum register the number $i/2^m$:

from qiskit.circuit.library import PiecewiseChebyshev 
from qiskit import *
from qiskit.algorithms import AmplitudeEstimation, EstimationProblem, AmplitudeEstimationResult
from qiskit.utils import QuantumInstance
import numpy as np
from qiskit.extensions.quantum_initializer.initializer import Initialize

def create_state_preparation_for_arcsin_sqrt(N, number):
    """N: input bits
       number: 2**N elements vectors whose only non zero component is the i-th one, and number represents |i> vector
    if len(number) != 2**N:
        raise ValueError("Number must be a 2**N elements vectors")
    func = lambda x: np.arcsin(np.sqrt(x / 2 ** N)) # The function to be implemented
    degree = 5 # The degree of Chebyshev polynomials. Use higher degree for better approximation
    breakpoints = [0, 2**N-1]
    pw_approx = PiecewiseChebyshev(func, degree, breakpoints, N)
    qr_state = QuantumRegister(N, "state")
    qr_output = QuantumRegister(1, "output")
    qr_aux = QuantumRegister(pw_approx.num_ancillas, "aux")

    qc = QuantumCircuit(N + 1 + pw_approx.num_ancillas)
    gate = Initialize(number).gates_to_uncompute().inverse()
    qc.append(gate, range(N))
    qc.append(pw_approx.to_instruction(), range(N+1+pw_approx.num_ancillas))
    qc2 = transpile(qc, basis_gates=["u3", "cx"])
    return qc2, [N]

N = 3
number = np.array([0] * (2**3))
number[3] = 1
state_preparation, objective_qubits = create_state_preparation_for_arcsin_sqrt(3, number)

backend = Aer.get_backend('qasm_simulator', shots=10000)
qinstance = QuantumInstance(backend, seed_simulator=2, seed_transpiler=2)
m = 4

problem = EstimationProblem(state_preparation=state_preparation, objective_qubits=objective_qubits)
ae = AmplitudeEstimation(num_eval_qubits=m, quantum_instance=qinstance)

result = ae.estimate(problem)
print('Grid-based estimate:', result.estimation)
print('Improved continuous estimate:', result.mle)

Once you have defined the AmplitudeEstimation object correctly you can append the corresponding circuit using

ae_circuit = ae.construct_circuit(estimation_problem=problem)

and use it inside your workflow. However, it seems cheaper to implements the quantum arithmetical gates implementing $|x\rangle_n |0...0\rangle_m \to |x\rangle_n|f(x)\rangle_m$ although it requires much more work.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.