Error with jnp array in multiplying it to a Hamiltonian

NickGut0711 · January 23, 2024, 3:54am

Hello!

I’m currently creating a QNode, cost, with the JAX interface.

def get_observables(N):
    
    observables = []

    # Coupling operators
    for i in range(N-1):
        observables.append(qml.PauliZ(i) @ qml.PauliZ(i+1))

    # Identity operator
    for i in range(N):
        observables.append(qml.Identity(i))

    return observables

def get_coeffs(params, N):

    coeffs = []

    # Coupling coeffs
    for i in range(N-1):
        coeffs.append((params[0])**2/params[1])

    # Constant coeffs
    for i in range(N):
        coeffs.append(params[1])

    return coeffs

def create_Hamiltonian(params):

    coeffs = get_coeffs(params, nqubits)
    obs = get_observables(nqubits)
    
    # H = qml.dot(coeffs, obs)
    H = qml.Hamiltonian(coeffs, obs)

    return H

def get_Sy(nqubits, a):
    S_0 = nqubits/2
    c = 0

    for i in range(nqubits):
        c += (1/(2*S_0))*qml.PauliY(wires=i)
    return a * c

def create_params(L, scale = 0.1):
    
    params = jnp.array([])

    print(type(params))

    for i in range(L):
        J = scale*np.random.uniform()
        O = 1.0
        theta = scale*np.random.uniform()

        # Convert the list to a NumPy array
        param_array = jnp.array([J, O, theta])
        
        params = jnp.append(params, param_array)

    print(type(params))
            
    return params

def U1(params):
    start_index = 0
    num_trotter_steps = 10

    for i in range(L):
        new_params = params[start_index:start_index + 3]
        H = create_Hamiltonian(new_params[0:2])
        qml.evolve(H, num_steps = num_trotter_steps)
        for j in range(nqubits):
            qml.RX(new_params[2], wires=j) # Change params to make sure that theta value changes for each L
        start_index += 3 # Put state_index in again

dev = qml.device("default.qubit.jax", wires=nqubits, shots=None)

# @jax.jit
# @qml.qnode(dev, interface='jax')
def circuit(params, a, phi):
    print(type(params))
    print(type(a))

    for i in range(nqubits): # Making the initial CSS
        qml.Hadamard(wires=i)

    U1(params)

    for z in range(nqubits): # Perturbation
        qml.RY(phi, wires = z)

    qml.adjoint(U1)(params)

    # expectation_values = [qml.expval(qml.PauliY(wires=i)) for i in range(nqubits)]

    c = get_Sy(nqubits, a)

    return qml.expval(c)

# @jax.jit
@qml.qnode(dev, interface='jax')
def cost(params, a, phi):
    circuit_output = circuit(params, a, phi)
    mse = np.mean((phi-circuit_output)**2)
    return mse

I set the parameters and am trying to calculate the gradients via jax.value_and_grad.

L = 1
params = create_params(L)
phi = jnp.array(0.001)
a = jnp.array(0.001)


val, grads = jax.value_and_grad(cost)(params, a, phi)

which gives me the error

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/.../.ipynb 
phi = jnp.array(0.001)
/.../.ipynb
a = jnp.array(0.001)
----> /.../.ipynb
val, grads = jax.value_and_grad(cost)(params, a, phi)

    [... skipping hidden 8 frame]

File ~/.../, in QNode.__call__(self, *args, **kwargs)
    967         kwargs[\"shots\"] = _get_device_shots(self._original_device)
    969 # construct the tape
--> 970 self.construct(args, kwargs)
    972 cache = self.execute_kwargs.get(\"cache\", False)
    973 using_custom_cache = (
    974     hasattr(cache, \"__getitem__\")
    975     and hasattr(cache, \"__setitem__\")
    976     and hasattr(cache, \"__delitem__\")
    977 )

File ~/.../ in QNode.construct(self, args, kwargs)
    853     self.interface = qml.math.get_interface(*args, *list(kwargs.values()))
    855 with qml.queuing.AnnotatedQueue() as q:
--> 856     self._qfunc_output = self.func(*args, **kwargs)
    858 self._tape = QuantumScript.from_queue(q, shots)
    860 params = self.tape.get_parameters(trainable_only=False)

/.../.ipynb Cell 9 line 1
  @qml.qnode(dev, interface='jax')
    def cost(params, a, phi):
--> circuit_output = circuit(params, a, phi)
    mse = np.mean((phi-circuit_output)**2)
    return mse

/.../.ipynb Cell 9 line 1
     qml.adjoint(U1)(params)
    # expectation_values = [qml.expval(qml.PauliY(wires=i)) for i in range(nqubits)]
-->  c = get_Sy(nqubits, a)
     return qml.expval(c)

/.../.ipynb Cell 9 line 4
for i in range(nqubits):
     c += (1/(2*S_0))*qml.PauliY(wires=i)
---> return a * c

TypeError: unsupported operand type(s) for *: 'ArrayImpl' and 'Hamiltonian'

This is what happens when I type qml.about().

Name: PennyLane
Version: 0.33.1
Summary: PennyLane is a Python quantum machine learning library by Xanadu Inc.
Home-page: https://github.com/PennyLaneAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: /.../python3.11/site-packages
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions
Required-by: PennyLane-Lightning

Platform info:           macOS-13.4.1-x86_64-i386-64bit
Python version:          3.11.5
Numpy version:           1.26.2
Scipy version:           1.11.4
Installed devices:
- default.gaussian (PennyLane-0.33.1)
- default.mixed (PennyLane-0.33.1)
- default.qubit (PennyLane-0.33.1)
- default.qubit.autograd (PennyLane-0.33.1)
- default.qubit.jax (PennyLane-0.33.1)
- default.qubit.legacy (PennyLane-0.33.1)
- default.qubit.tf (PennyLane-0.33.1)
- default.qubit.torch (PennyLane-0.33.1)
- default.qutrit (PennyLane-0.33.1)
- null.qubit (PennyLane-0.33.1)
- lightning.qubit (PennyLane-Lightning-0.33.1)

I’m unsure what it is that is causing this error since I’m doing similar things with params as I am with a inside the circuit QNode.

isaacdevlugt · January 24, 2024, 2:31pm

Hey @NickGut0711,

First off, default.qubit.jax uses the old device API. If you just use default.qubit (new device API) it will automatically see what interface to use

When multiplying jax scalars by Hamiltonians, you should use qml.s_prod: qml.s_prod — PennyLane 0.34.0 documentation. That should take care of your error! Let me know if that’s the case

NickGut0711 · January 24, 2024, 3:38pm

Hi @isaacdevlugt! Thank you for your response. I actually used qml.s_prod when I was trying to figure this out on my own, but I was using the old device API you mentioned. I kept getting an error like:

   1323     eigvals = self._asarray(
   1324         observable.eigvals()
   1325         if not isinstance(observable, MeasurementValue)
   1326         else np.arange(2 ** len(observable.wires)),
   1327         dtype=self.R_DTYPE,
   1328     )
   1329 except qml.operation.EigvalsUndefinedError as e:
-> 1330     raise qml.operation.EigvalsUndefinedError(
   1331         f"Cannot compute analytic expectations of {observable.name}."
   1332     ) from e
   1334 prob = self.probability(wires=observable.wires)
   1335 # In case of broadcasting, `prob` has two axes and this is a matrix-vector product

EigvalsUndefinedError: Cannot compute analytic expectations of SProd.

Now that I updated the device API, and implemented qml.s_prod like I do below:

def get_observables(N):
    
    observables = []

    # Coupling operators
    for i in range(N-1):
        observables.append(qml.PauliZ(i) @ qml.PauliZ(i+1))

    # Identity operator
    for i in range(N):
        observables.append(qml.Identity(i))

    return observables

def get_coeffs(params, N):

    coeffs = []

    # Coupling coeffs
    for i in range(N-1):
        coeffs.append((params[0])**2/params[1])

    # Constant coeffs
    for i in range(N):
        coeffs.append(params[1])

    return coeffs

def create_Hamiltonian(params):

    coeffs = get_coeffs(params, nqubits)
    obs = get_observables(nqubits)
    
    # H = qml.dot(coeffs, obs)
    H = qml.Hamiltonian(coeffs, obs)

    return H

def get_Sy(nqubits, a):
    S_0 = nqubits/2
    c = 0

    for i in range(nqubits):
        c += (1/(2*S_0))*qml.PauliY(wires=i)
    
    return qml.s_prod(a, c)

# def get_Sy(nqubits):
#     S_0 = nqubits / 2
#     pauli_terms = [qml.PauliY(wires=i) for i in range(nqubits)]
#     hamiltonian = qml.operation.Tensor(*pauli_terms) / (2 * S_0)
#     return hamiltonian

def create_params(L, scale = 0.1):
    
    params = jnp.array([])

    print('Params in create_params before:' + str(type(params)))

    for i in range(L):
        J = scale*np.random.uniform()
        O = 1.0
        theta = scale*np.random.uniform()

        # Convert the list to a NumPy array
        param_array = jnp.array([J, O, theta])
        
        params = jnp.append(params, param_array)

    print('Params in create_params after:' + str(type(params)))
            
    return params

def U1(params):
    start_index = 0
    num_trotter_steps = 10

    for i in range(L):
        new_params = params[start_index:start_index + 3]
        H = create_Hamiltonian(new_params[0:2])
        qml.evolve(H, num_steps = num_trotter_steps)
        for j in range(nqubits):
            qml.RX(new_params[2], wires=j) # Change params to make sure that theta value changes for each L
        start_index += 3 # Put state_index in again

dev = qml.device("default.qubit", wires=nqubits, shots=None)

# @jax.jit
@qml.qnode(dev, interface='jax')
def circuit(params, a, phi):

    for i in range(nqubits): # Making the initial CSS
        qml.Hadamard(wires=i)

    U1(params)

    for z in range(nqubits): # Perturbation
        qml.RY(phi, wires = z)

    qml.adjoint(U1)(params)

    # expectation_values = [qml.expval(qml.PauliY(wires=i)) for i in range(nqubits)]

    c = get_Sy(nqubits, a)

    return qml.expval(c)

def mse(expected, predictions):
    loss = 0
    for l, p in zip(expected, predictions):
        loss = loss + (l - p)**2

    loss = loss / len(expected)
    
    return loss

def cost(params, a, expected):
    preds = [circuit(params, a, i) for i in range(len(expected))]
    return mse(expected, preds)

L = 1
params = create_params(L)
phi = jnp.array(0.001)
a = jnp.array(0.001)

circuit(params, a, phi)

I get a strange error that I have no idea how to fix!

NotImplementedError                       Traceback (most recent call last)
/.../.ipynb Cell 6 line 1
----> 1 circuit(params, a, phi)

File ~/.../pennylane/qnode.py:1027, in QNode.__call__(self, *args, **kwargs)
   1022         full_transform_program._set_all_argnums(
   1023             self, args, kwargs, argnums
   1024         )  # pylint: disable=protected-access
   1026 # pylint: disable=unexpected-keyword-arg
-> 1027 res = qml.execute(
   1028     (self._tape,),
   1029     device=self.device,
   1030     gradient_fn=self.gradient_fn,
   1031     interface=self.interface,
   1032     transform_program=full_transform_program,
   1033     config=config,
   1034     gradient_kwargs=self.gradient_kwargs,
   1035     override_shots=override_shots,
   1036     **self.execute_kwargs,
   1037 )
   1039 res = res[0]
   1041 # convert result to the interface in case the qfunc has no parameters

File ~/.../pennylane/interfaces/execution.py:616, in execute(tapes, device, gradient_fn, interface, transform_program, config, grad_on_execution, gradient_kwargs, cache, cachesize, max_diff, override_shots, expand_fn, max_expansion, device_batch_transform)
    614 # Exiting early if we do not need to deal with an interface boundary
    615 if no_interface_boundary_required:
--> 616     results = inner_execute(tapes)
    617     return post_processing(results)
    619 _grad_on_execution = False

File ~/.../pennylane/interfaces/execution.py:249, in _make_inner_execute.<locals>.inner_execute(tapes, **_)
    247 if numpy_only:
    248     tapes = tuple(qml.transforms.convert_to_numpy_parameters(t) for t in tapes)
--> 249 return cached_device_execution(tapes)

File ~/.../pennylane/interfaces/execution.py:371, in cache_execute.<locals>.wrapper(tapes, **kwargs)
    366         return (res, []) if return_tuple else res
    368 else:
    369     # execute all unique tapes that do not exist in the cache
    370     # convert to list as new device interface returns a tuple
--> 371     res = list(fn(tuple(execution_tapes.values()), **kwargs))
    373 final_res = []
    375 for i, tape in enumerate(tapes):

File ~/.../pennylane/devices/default_qubit.py:474, in DefaultQubit.execute(self, circuits, execution_config)
    468 interface = (
    469     execution_config.interface
    470     if execution_config.gradient_method in {"backprop", None}
    471     else None
    472 )
    473 if max_workers is None:
--> 474     results = tuple(
    475         simulate(
    476             c,
    477             rng=self._rng,
    478             prng_key=self._prng_key,
    479             debugger=self._debugger,
    480             interface=interface,
    481         )
    482         for c in circuits
    483     )
    484 else:
    485     vanilla_circuits = [convert_to_numpy_parameters(c) for c in circuits]

File ~/.../pennylane/devices/default_qubit.py:475, in <genexpr>(.0)
    468 interface = (
    469     execution_config.interface
    470     if execution_config.gradient_method in {"backprop", None}
    471     else None
    472 )
    473 if max_workers is None:
    474     results = tuple(
--> 475         simulate(
    476             c,
    477             rng=self._rng,
    478             prng_key=self._prng_key,
    479             debugger=self._debugger,
    480             interface=interface,
    481         )
    482         for c in circuits
    483     )
    484 else:
    485     vanilla_circuits = [convert_to_numpy_parameters(c) for c in circuits]

File ~/.../pennylane/devices/qubit/simulate.py:270, in simulate(circuit, rng, prng_key, debugger, interface)
    241 """Simulate a single quantum script.
    242 
    243 This is an internal function that will be called by the successor to ``default.qubit``.
   (...)
    267 
    268 """
    269 state, is_state_batched = get_final_state(circuit, debugger=debugger, interface=interface)
--> 270 return measure_final_state(circuit, state, is_state_batched, rng=rng, prng_key=prng_key)

File ~/.../pennylane/devices/qubit/simulate.py:211, in measure_final_state(circuit, state, is_state_batched, rng, prng_key)
    207 if not circuit.shots:
    208     # analytic case
    210     if len(circuit.measurements) == 1:
--> 211         return measure(circuit.measurements[0], state, is_state_batched=is_state_batched)
    213     return tuple(
    214         measure(mp, state, is_state_batched=is_state_batched) for mp in circuit.measurements
    215     )
    217 # finite-shot case

File ~/.../pennylane/devices/qubit/measure.py:197, in measure(measurementprocess, state, is_state_batched)
    184 def measure(
    185     measurementprocess: MeasurementProcess, state: TensorLike, is_state_batched: bool = False
    186 ) -> TensorLike:
    187     """Apply a measurement process to a state.
    188 
    189     Args:
   (...)
    195         Tensorlike: the result of the measurement
    196     """
--> 197     return get_measurement_function(measurementprocess, state)(
    198         measurementprocess, state, is_state_batched
    199     )

File ~/.../pennylane/devices/qubit/measure.py:181, in get_measurement_function(measurementprocess, state)
    178     if measurementprocess.obs is None or measurementprocess.obs.has_diagonalizing_gates:
    179         return state_diagonalizing_gates
--> 181 raise NotImplementedError

NotImplementedError:

isaacdevlugt · January 26, 2024, 11:59pm

Interesting . I’m a little stumped here myself. Can you submit an issue to our github repo?

Topic		Replies	Views
How do I encode the MSE in a certain cost function? PennyLane Help	1	184	January 24, 2024
QubitUnitary and jax.jit PennyLane Help	2	353	February 28, 2022
Vectorizing pennylane neural networks PennyLane Help	1	14	August 12, 2024
From `qml.Hamiltonian` to numpy representation PennyLane Help	12	2020	March 20, 2024
Qml.counts is not working with jax.jit PennyLane Help	6	608	March 6, 2023

Error with jnp array in multiplying it to a Hamiltonian

Related topics