Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Amend PCPhase operator to support RZ convention #3926

Merged
merged 7 commits into from
Apr 3, 2023
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 28 additions & 16 deletions pennylane/ops/qubit/parametric_ops_multi_qubit.py
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
Expand Up @@ -558,13 +558,13 @@ class PCPhase(Operation):

basis = "Z"
grad_method = "A"
parameter_frequencies = [(1,)]
parameter_frequencies = [(2,)]
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved

def generator(self):
r"""The hermitian matrix, which when exponentiated and scaled (i :math:`\phi`), produces
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
the PCPhase unitary."""
dim, shape = self.hyperparameters["dimension"]
mat = np.diag([1 if index < dim else 0 for index in range(shape)])
mat = np.diag([1 if index < dim else -1 for index in range(shape)])
return qml.Hermitian(mat, wires=self.wires)

def __init__(self, phi, dim, wires, do_queue=True, id=None):
Expand All @@ -587,21 +587,21 @@ def compute_matrix(*params, **hyperparams):

if qml.math.get_interface(phi) == "tensorflow":
p = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
ones = qml.math.ones_like(p)
minus_p = qml.math.exp(-1j * qml.math.cast_like(phi, 1j))
zeros = qml.math.zeros_like(p)

columns = []
for i in range(t):
columns.append(
[p if j == i else zeros for j in range(t)]
if i < d
else [ones if j == i else zeros for j in range(t)]
else [minus_p if j == i else zeros for j in range(t)]
)
r = qml.math.stack(columns, like="tensorflow", axis=-2)
return r

arg = 1j * phi
prefactors = qml.math.array([1 if index < d else 0 for index in range(t)], like=phi)
prefactors = qml.math.array([1 if index < d else -1 for index in range(t)], like=phi)

if qml.math.ndim(arg) == 0:
return qml.math.diag(qml.math.exp(arg * prefactors))
Expand All @@ -617,12 +617,11 @@ def compute_eigvals(*params, **hyperparams):

if qml.math.get_interface(phi) == "tensorflow":
phase = qml.math.exp(1j * qml.math.cast_like(phi, 1j))
return stack_last(
[phase if index < d else qml.math.ones_like(phase) for index in range(t)]
)
minus_phase = qml.math.exp(-1j * qml.math.cast_like(phi, 1j))
return stack_last([phase if index < d else minus_phase for index in range(t)])

arg = 1j * phi
prefactors = qml.math.array([1 if index < d else 0 for index in range(t)], like=phi)
prefactors = qml.math.array([1 if index < d else -1 for index in range(t)], like=phi)

if qml.math.ndim(phi) == 0:
product = arg * prefactors
Expand Down Expand Up @@ -654,23 +653,36 @@ def compute_decomposition(*params, wires=None, **hyperparams):
phi = params[0]
k, n = hyperparams["dimension"]

def _get_base_ops(theta, wire):
return PauliX(wire) @ PhaseShift(theta, wire) @ PauliX(wire), PhaseShift(theta, wire)

def _get_op_from_binary_rep(binary_rep, theta, wires):
if len(binary_rep) == 1:
op = _get_base_ops(theta, wire=wires[0])[int(binary_rep)]
op = (
PhaseShift(theta, wires[0])
if int(binary_rep)
else PauliX(wires[0]) @ PhaseShift(theta, wires[0]) @ PauliX(wires[0])
)
else:
base_op = _get_base_ops(theta, wire=wires[-1])[int(binary_rep[-1])]
base_op = (
PhaseShift(theta, wires[-1])
if int(binary_rep[-1])
else PauliX(wires[-1]) @ PhaseShift(theta, wires[-1]) @ PauliX(wires[-1])
)
op = qml.ctrl(
base_op, control=wires[:-1], control_values=[int(i) for i in binary_rep[:-1]]
)
return op

n_log2 = int(np.log2(n))
binary_reps = [bin(_k)[2:].zfill(n_log2) for _k in range(k)]
positive_binary_reps = [bin(_k)[2:].zfill(n_log2) for _k in range(k)]
negative_binary_reps = [bin(_k)[2:].zfill(n_log2) for _k in range(k, n)]

positive_ops = [
_get_op_from_binary_rep(br, phi, wires=wires) for br in positive_binary_reps
]
negative_ops = [
_get_op_from_binary_rep(br, -1 * phi, wires=wires) for br in negative_binary_reps
]

return [_get_op_from_binary_rep(br, phi, wires=wires) for br in binary_reps]
return positive_ops + negative_ops

def adjoint(self):
"""Computes the adjoint of the operator."""
Expand Down
2 changes: 1 addition & 1 deletion tests/ops/qubit/test_attributes.py
Original file line number Diff line number Diff line change
Expand Up @@ -465,7 +465,7 @@ def test_pcphase(self):
mat2 = op.compute_matrix(*op.parameters, **op.hyperparameters)

mats = [
np.diag([np.exp(1j * phi) if i < dim else 1.0 for i in range(size)])
np.diag([np.exp(1j * phi) if i < dim else np.exp(-1j * phi) for i in range(size)])
for phi in broadcasted_phi
]
expected_mat = np.array(mats)
Expand Down
48 changes: 29 additions & 19 deletions tests/ops/qubit/test_parametric_ops.py
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,13 @@ def circuit(phase, dim):
return qml.state()

def _get_expected_state(phase, dim, size):
return np.array([np.exp(1j * phase) if i < dim else 1.0 for i in range(size)]) * 1 / 2
return (
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
np.array(
[np.exp(1j * phase) if i < dim else np.exp(-1j * phase) for i in range(size)]
)
* 1
/ 2
)

assert np.allclose(circuit(theta, d), _get_expected_state(theta, d, 4))

Expand Down Expand Up @@ -996,7 +1002,7 @@ def test_pcphase(self, phi, dim, wires):
mat2 = op.compute_matrix(*op.parameters, **op.hyperparameters)

expected_mat = np.diag(
[np.exp(1j * phi) if i < dim else 1.0 for i in range(2**num_wires)]
[np.exp(1j * phi) if i < dim else np.exp(-1j * phi) for i in range(2**num_wires)]
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
)
assert np.allclose(mat1, expected_mat)
assert np.allclose(mat2, expected_mat)
Expand All @@ -1017,7 +1023,9 @@ def test_pcphase_tf(self, phi, dim, wires):
mat2 = op.compute_matrix(*op.parameters, **op.hyperparameters)

expected_mat = tf.Variable(
np.diag([np.exp(1j * phi) if i < dim else 1.0 for i in range(2**num_wires)])
np.diag(
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
[np.exp(1j * phi) if i < dim else np.exp(-1j * phi) for i in range(2**num_wires)]
)
)

assert np.allclose(mat1, expected_mat)
Expand All @@ -1038,7 +1046,9 @@ def test_pcphase_torch(self, phi, dim, wires):
mat2 = op.compute_matrix(*op.parameters, **op.hyperparameters)

expected_mat = torch.tensor(
np.diag([np.exp(1j * phi) if i < dim else 1.0 for i in range(2**num_wires)])
np.diag(
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
[np.exp(1j * phi) if i < dim else np.exp(-1j * phi) for i in range(2**num_wires)]
)
)

assert np.allclose(mat1, expected_mat)
Expand All @@ -1061,7 +1071,9 @@ def test_pcphase_jax(self, phi, dim, wires):
mat2 = op.compute_matrix(*op.parameters, **op.hyperparameters)

expected_mat = jnp.diag(
jnp.array([jnp.exp(1j * phi) if i < dim else 1.0 for i in range(2**num_wires)])
jnp.array(
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
[jnp.exp(1j * phi) if i < dim else np.exp(-1j * phi) for i in range(2**num_wires)]
)
)

assert np.allclose(mat1, expected_mat)
Expand All @@ -1080,7 +1092,7 @@ def test_pcphase_broadcasted(self):
mat2 = op.compute_matrix(*op.parameters, **op.hyperparameters)

mats = [
np.diag([np.exp(1j * phi) if i < dim else 1.0 for i in range(size)])
np.diag([np.exp(1j * phi) if i < dim else np.exp(-1j * phi) for i in range(size)])
for phi in broadcasted_phi
]
expected_mat = np.array(mats)
Expand Down Expand Up @@ -2105,13 +2117,17 @@ def test_pcphase_eigvals(self):
# test arbitrary phase shift
phi = 0.5432
op = qml.PCPhase(phi, dim=2, wires=[0, 1])
expected = np.array([np.exp(1j * phi), np.exp(1j * phi), 1.0, 1.0])
expected = np.array(
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
[np.exp(1j * phi), np.exp(1j * phi), np.exp(-1j * phi), np.exp(-1j * phi)]
)
assert np.allclose(op.eigvals(), expected)

# test arbitrary broadcasted phase shift
phi = np.array([0.5, 0.4, 0.3])
op = qml.PCPhase(phi, dim=2, wires=[0, 1])
expected = np.array([[np.exp(1j * p), np.exp(1j * p), 1.0, 1.0] for p in phi])
expected = np.array(
Jaybsoni marked this conversation as resolved.
Show resolved Hide resolved
[[np.exp(1j * p), np.exp(1j * p), np.exp(-1j * p), np.exp(-1j * p)] for p in phi]
)
assert np.allclose(op.eigvals(), expected)

def test_crz_eigvals(self, tol):
Expand Down Expand Up @@ -2805,7 +2821,7 @@ def test_pcphase_grad(self, dev_name, diff_method, phi):
pytest.skip("PCPhase does not support adjoint diff")

dev = qml.device(dev_name, wires=[0, 1])
expected_grad = (1j / 2) * (npp.exp(1j * phi) - npp.exp(-1j * phi)) # computed by hand
expected_grad = -4 * npp.cos(phi) * npp.sin(phi) # computed by hand

@qml.qnode(dev, diff_method=diff_method)
def circ(phi):
Expand All @@ -2832,9 +2848,7 @@ def test_pcphase_grad_tf(self, dev_name, diff_method, phi):
import tensorflow as tf

dev = qml.device(dev_name, wires=[0, 1])
expected_grad = tf.Variable(
(1j / 2) * (npp.exp(1j * phi) - npp.exp(-1j * phi))
) # computed by hand
expected_grad = tf.Variable(-4 * npp.cos(phi) * npp.sin(phi)) # computed by hand
phi = tf.Variable(phi)

@qml.qnode(dev, diff_method=diff_method)
Expand Down Expand Up @@ -2864,9 +2878,7 @@ def test_pcphase_grad_torch(self, dev_name, diff_method, phi):
import torch

dev = qml.device(dev_name, wires=[0, 1])
expected_grad = torch.tensor(
(1j / 2) * (npp.exp(1j * phi) - npp.exp(-1j * phi))
) # computed by hand
expected_grad = torch.tensor(-4 * npp.cos(phi) * npp.sin(phi)) # computed by hand
phi = torch.tensor(phi, requires_grad=True)

@qml.qnode(dev, diff_method=diff_method)
Expand Down Expand Up @@ -2895,9 +2907,7 @@ def test_pcphase_grad_jax(self, dev_name, diff_method, phi):
import jax.numpy as jnp

dev = qml.device(dev_name, wires=[0, 1])
expected_grad = jnp.array(
(1j / 2) * (npp.exp(1j * phi) - npp.exp(-1j * phi))
) # computed by hand
expected_grad = jnp.array(-4 * npp.cos(phi) * npp.sin(phi)) # computed by hand
phi = jnp.array(phi)

@qml.qnode(dev, diff_method=diff_method)
Expand Down Expand Up @@ -2939,7 +2949,7 @@ def test_pcphase_generator(self):
phi = 1.23
op = qml.PCPhase(phi, dim=2, wires=[0, 1])

expected_mat = np.diag([1.0, 1.0, 0.0, 0.0])
expected_mat = np.diag([1.0, 1.0, -1.0, -1.0])

gen, coeff = qml.generator(op)
assert np.allclose(qml.matrix(gen), expected_mat)
Expand Down