Arbitrary State preparation


I’m currently working on the Variational Quantum Linear Solver tutorial, and I would like to generalize it to any vector b and matrix A.

Is there anyway to do that with pennylane?

Best regards

Hi @Ramy! You can definitely do that.

You would need to redefine the function U_b and CA defined in the Circuits of the quantum linear problem section.

You will probably also need to edit the variational_block function. The StronglyEntanglingLayers is a good place to start.

And finally, you will need to change the hyperparameters, with the new number of qubits.

Let me know if this works :slight_smile:

Hi @christina,
Thanks for your response!

I know that these functions need to be edited (U_b and CA) but I’m still new to pennylane so I don’t know how to do that. Is there a function similar to initialize function in qiskit which prepares any arbitrary state?

I’d recommend taking a look at Mottonen state preperation and the control transform.

1 Like

Thank you!
I will check that

1 Like

Hi @christina again,

For the state preparation, I think it is clear now, but I’m still confused with the matrix part. Suppose I have a hermitian matrix and I would like to convert it into operations, Do you have any idea on how to do it?


Hi @Ramy!

Are you looking for the qml.utils.decompose_hamiltonian function?

This function takes an Hermitian matrix, and returns a tuple of coefficients and Pauli operators:

>>> A = np.array( ... [[-2, -2+1j, -2, -2], [-2-1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]])
>>> coeffs, obs_list = decompose_hamiltonian(A)
>>> coeffs [-1.0, -1.5, -0.5, -1.0, -1.5, -1.0, -0.5, 1.0, -0.5, -0.5]

We can use the output coefficients and tensor Pauli terms to construct a Hamiltonian :

>>> H = qml.Hamiltonian(coeffs, obs_list)
>>> print(H)
(-1.0) [I0 I1]
 + (-1.5) [X1]
 + (-0.5) [Y1]
 + (-1.0) [Z1]
 + (-1.5) [X0]
 + (-1.0) [X0 X1]
 + (-0.5) [X0 Z1]
 + (1.0) [Y0 Y1]
 + (-0.5) [Z0 X1]
 + (-0.5) [Z0 Y1]
1 Like

Thank you for your response, I think this is what I need.
just to make sure, it should work with qml.ctrl right?

Hi @Ramy, you might have to provide more details in terms of how you would like to use it with qml.ctrl :slight_smile:

Hi @josh, I’m trying to generalize Variational Quantum Linear for any Hermitian matrix A and arbitrary state b.

Hi @Ramy,

Thanks for the questions! qml.ctrl is mainly meant to be used with quantum operations, rather than an observable or Hamiltonian. How could qml.ctrl fit into your use case of generalizing the Variational Quantum Linear solver?

1 Like

Hi @antalszava, thanks for the clarification!
Actually, I’m not sure if it can be used or not.
My original question was how to generalize the Variational Quantum Linear solver for any hermition matrix?

Hi @Ramy,

The previously mentioned components are useful in achieving that.

  1. For changing U_b, the proposed qml.MottonenStatePreparation can be used.
  2. For changing CA, first let’s have a look at how we can decompose an Hermitian matrix A. We can represent A as A = \sum_{I}C_I\hat{P}_I

where C_I are numerical coefficients, and P_I are Pauli
“words”, products of Pauli operators of different qubits \hat{P}_I=\prod_{i=1}^N\hat{\sigma}_i^{(I)}


\hat{\sigma}_i^{(I)} is one of the \hat{x}, \hat{y}, \hat{z} Pauli operators or identity \hat{e} for the i^{th} qubit.

(As described in the introduction of this reference for qubit Hamiltonians, i.e., Hermitian operators).

To get this decomposition, the decompose_hamiltonian function mentioned by Josh can be used as basis. The main difference, from how this function is being used usually, is that we’d create a Hamiltonian class using qml.Hamiltonian. A qml.Hamiltonian object is usually used with qml.ExpvalCost to return the expectation value of the Pauli words.

In the tutorial, however, we are applying the unitaries (i.e., the Pauli words from our decomposition) as quantum operations.

Therefore, we can use a custom definition of the decompose_hamiltonian function:

def decompose_hamiltonian_unitary_ops(H):
    """Auxiliary function similar to pennylane.utils.decompose_hamiltonian.
    This function returns a nested list of 
    n = int(np.log2(len(H)))
    N = 2 ** n

    if H.shape != (N, N):
        raise ValueError(
            "The Hamiltonian should have shape (2**n, 2**n), for any qubit number n>=1"

    if not np.allclose(H, H.conj().T):
        raise ValueError("The Hamiltonian is not Hermitian")

    paulis = [qml.Identity, qml.PauliX, qml.PauliY, qml.PauliZ]
    ops = []
    coeffs = []

    for term in itertools.product(paulis, repeat=n):
        matrices = [i._matrix() for i in term]
        coeff = np.trace(functools.reduce(np.kron, matrices) @ H) / N
        coeff = np.real_if_close(coeff).item()

        if not np.allclose(coeff, 0):

            if not all(t is qml.Identity for t in term) and True:
                ops.append([t(i) for i, t in enumerate(term) if t is not qml.Identity])

    return coeffs, ops

This definition is very similar to the original implementation. The key difference is, that with the original version we are creating the tensor products of Paulis, whereas with this version, we accumulate them in lists that are contained in ops.

Once we have this function, we can decompose arbitrary Hermitian matrices:

H = np.array([[-2, -2+1j, -2, -2], [-2-1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]])

coeffs, ops_list = decompose_hamiltonian_unitary_ops(H)

At this point, we just have to define a quantum function (e.g., called applying_unitaries), that will apply the unitaries that we’ve gathered in ops_list. Under the hood, PennyLane keeps track of operations in a quantum function using a concept called queuing. Each operation in the circuit is being queued when an operation is created when calling the quantum function.

For example, the qml.RY and the qml.PauliZ operators are queued in the following example (this example is not required for the solution):

def circuit():
    qml.RY(0.3, wires=[0])
    return qml.expval(qml.PauliZ(0))

In our case, however, we have already created our operations, so we have to explicitly queue them:

def applying_unitaries():
    """Quantum function that applies the unitaries building up a Hamiltonian"""
    for ops in ops_list:
        for o in ops:

Once we have this quantum function, all is ready to make a controlled version of it using qml.ctrl:

# Creating a controlled version for applying the unitaries
# The second argument are the control wire(s) to use
ctrl_ops = qml.ctrl(applying_unitaries, [2])

ctrl_ops will return a quantum function that can be used in any QNode. We can check that we have the correct circuit with controlled operations by comparing the original circuit and the new circuit:

dev = qml.device('default.qubit', 3)

def applying_unitaries_with_expval():
    return qml.expval(qml.PauliZ(0))

def ctrl_circuit_with_expval():
    return qml.expval(qml.PauliZ(0))


 0: ──X──X──X──Y──Z──Z────────┤ ⟨Z⟩ 
 1: ──X──Y──Z──X──Z──Y──X──Y──┤     

 0: ───────────────╭X──╭X──────╭X──────╭CY───────╭Z──────╭Z───────┤ ⟨Z⟩ 
 1: ──╭X──╭CY──╭Z──│───│───╭X──│───╭Z──│────╭CY──│───╭X──│───╭CY──┤     
 2: ──╰C──╰CY──╰C──╰C──╰C──╰C──╰C──╰C──╰CY──╰CY──╰C──╰C──╰C──╰CY──┤     

To adjust this code to be written as CA in the tutorial, we will likely need to index into the list of operations.

Hope this helps! :slightly_smiling_face:

1 Like

Thank you very much for the detailed explanation @antalszava!
It will take me some time to go through it but I think it is clear now.
I appreciate it :blush:

1 Like

Hi @josh,
It’s not mentioned in VQLS implementation that A matrix needs to be hermitian as such but here, we can’t decompose it if it’s not hermitian. Do we need to take any other method to solve it when A isn’t hermitian or the algorithm needs a Hermitian matrix only? Please could you clarify?

I think there is some confusion between unitary and hermitian in this thread.

All “gates” in quantum computing are “unitary”.

Any “observable” or measurement made must be “hermitian”.

For the VQLS tutorial, both U_b and CA are unitary.

Josh and Antal provided excellent ways to decompose a hermitian matrix as per your question, but U_b and CA do not necessarily fit into that category.

The goal of U_b is to prepare a specific state. Arbitrary state preparation can be done with:

The goal of CA is to apply a specific controlled unitary. Given a unitary matrix, you can apply CA with qml.ControlledQubitUnitary

For smaller unitary matrix sizes, PennyLane can decompose a unitary into standard base gates like rotations using a variety of formulas.

Ya got the point, but for CA to be broken into Unitaries, the method I know is the same one just leaving the check which you have placed i.e., by calculating trace only. Is my understanding correct or are there other methods to decompose CA into linear combinations of Pauli Operators?

If a unitary is a linear combination of other operators, involving a sum, there’s no easy way to just apply it as the terms in the series. Instead, you have to calculate the effective matrix and then apply that.

Each gate performs a matrix multiplication on the state:

U_3 U_2 U_1 | 0 \rangle

Summing gates together like:

c_1 U_1 + c_2 U_2 = U_4

only creates a new matrix that we can instead apply.

To actually apply a sum of operators in a quantum circuit, not just with pen and paper, we would have to make a copy of the quantum state:

\left( c_1 U_1 + c_2 U_2 \right) | 0 \rangle = \left( c_1 U_1 |0\rangle \right) + \left( c_2 U_2 |0\rangle \right)

While we can do this on pen and paper, or on some simulator, we can’t do it on a quantum computer. And because we can’t do that on an actual quantum computer, we don’t support it on our simulators either.

I got your point, but the way of decomposition is as mentioned in this StackOverflow answer right? I see the same in your hermitian decomposition function. So I want to just confirm that the above method is the only one to find decompositions? Please, if you can check it and revert to me, it would be really helpful.

There is no only method to decompose matrices into a linear combination. Some will work all the time, but not be very efficient. Others work on a subset of cases.

For decomposing an arbitrary matrix into gates that we can apply in a circuit, you may want to investigate qml.transforms.unitary_to_rot — PennyLane 0.19.1 documentation.