# Error with jnp array in multiplying it to a Hamiltonian

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

U1(params)

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

# 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)

``````

which gives me the error

``````---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/.../.ipynb
phi = jnp.array(0.001)
/.../.ipynb
a = jnp.array(0.001)
----> /.../.ipynb

[... 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
# 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:
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.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.

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

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

U1(params)

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

# 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,
1031     interface=self.interface,
1032     transform_program=full_transform_program,
1033     config=config,
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)

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:
``````

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