/
operator.py
535 lines (457 loc) · 19.9 KB
/
operator.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
# -*- coding: utf-8 -*-
# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2019.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
Matrix Operator class.
"""
import copy
import re
from numbers import Number
import numpy as np
from qiskit.circuit.quantumcircuit import QuantumCircuit
from qiskit.circuit.instruction import Instruction
from qiskit.circuit.library.standard_gates import IGate, XGate, YGate, ZGate, HGate, SGate, TGate
from qiskit.exceptions import QiskitError
from qiskit.quantum_info.operators.predicates import is_unitary_matrix, matrix_equal
from qiskit.quantum_info.operators.base_operator import BaseOperator
class Operator(BaseOperator):
r"""Matrix operator class
This represents a matrix operator :math:`M` that will
:meth:`~Statevector.evolve` a :class:`Statevector` :math:`|\psi\rangle`
by matrix-vector multiplication
.. math::
|\psi\rangle \mapsto M|\psi\rangle,
and will :meth:`~DensityMatrix.evolve` a :class:`DensityMatrix` :math:`\rho`
by left and right multiplication
.. math::
\rho \mapsto M \rho M^\dagger.
"""
def __init__(self, data, input_dims=None, output_dims=None):
"""Initialize an operator object.
Args:
data (QuantumCircuit or
Instruction or
BaseOperator or
matrix): data to initialize operator.
input_dims (tuple): the input subsystem dimensions.
[Default: None]
output_dims (tuple): the output subsystem dimensions.
[Default: None]
Raises:
QiskitError: if input data cannot be initialized as an operator.
Additional Information:
If the input or output dimensions are None, they will be
automatically determined from the input data. If the input data is
a Numpy array of shape (2**N, 2**N) qubit systems will be used. If
the input operator is not an N-qubit operator, it will assign a
single subsystem with dimension specified by the shape of the input.
"""
if isinstance(data, (list, np.ndarray)):
# Default initialization from list or numpy array matrix
self._data = np.asarray(data, dtype=complex)
elif isinstance(data, (QuantumCircuit, Instruction)):
# If the input is a Terra QuantumCircuit or Instruction we
# perform a simulation to construct the unitary operator.
# This will only work if the circuit or instruction can be
# defined in terms of unitary gate instructions which have a
# 'to_matrix' method defined. Any other instructions such as
# conditional gates, measure, or reset will cause an
# exception to be raised.
self._data = self._init_instruction(data).data
elif hasattr(data, 'to_operator'):
# If the data object has a 'to_operator' attribute this is given
# higher preference than the 'to_matrix' method for initializing
# an Operator object.
data = data.to_operator()
self._data = data.data
if input_dims is None:
input_dims = data._input_dims
if output_dims is None:
output_dims = data._output_dims
elif hasattr(data, 'to_matrix'):
# If no 'to_operator' attribute exists we next look for a
# 'to_matrix' attribute to a matrix that will be cast into
# a complex numpy matrix.
self._array = np.asarray(data.to_matrix(), dtype=complex)
else:
raise QiskitError("Invalid input data format for Operator")
# Determine input and output dimensions
dout, din = self._data.shape
output_dims = self._automatic_dims(output_dims, dout)
input_dims = self._automatic_dims(input_dims, din)
super().__init__(input_dims, output_dims)
def __repr__(self):
prefix = 'Operator('
pad = len(prefix) * ' '
return '{}{},\n{}input_dims={}, output_dims={})'.format(
prefix, np.array2string(
self.data, separator=', ', prefix=prefix),
pad, self._input_dims, self._output_dims)
def __eq__(self, other):
"""Test if two Operators are equal."""
if not super().__eq__(other):
return False
return np.allclose(
self.data, other.data, rtol=self.rtol, atol=self.atol)
@property
def data(self):
"""Return data."""
return self._data
@classmethod
def from_label(cls, label):
"""Return a tensor product of single-qubit operators.
Args:
label (string): single-qubit operator string.
Returns:
Operator: The N-qubit operator.
Raises:
QiskitError: if the label contains invalid characters, or the
length of the label is larger than an explicitly
specified num_qubits.
Additional Information:
The labels correspond to the single-qubit matrices:
'I': [[1, 0], [0, 1]]
'X': [[0, 1], [1, 0]]
'Y': [[0, -1j], [1j, 0]]
'Z': [[1, 0], [0, -1]]
'H': [[1, 1], [1, -1]] / sqrt(2)
'S': [[1, 0], [0 , 1j]]
'T': [[1, 0], [0, (1+1j) / sqrt(2)]]
'0': [[1, 0], [0, 0]]
'1': [[0, 0], [0, 1]]
'+': [[0.5, 0.5], [0.5 , 0.5]]
'-': [[0.5, -0.5], [-0.5 , 0.5]]
'r': [[0.5, -0.5j], [0.5j , 0.5]]
'l': [[0.5, 0.5j], [-0.5j , 0.5]]
"""
# Check label is valid
label_mats = {
'I': IGate().to_matrix(),
'X': XGate().to_matrix(),
'Y': YGate().to_matrix(),
'Z': ZGate().to_matrix(),
'H': HGate().to_matrix(),
'S': SGate().to_matrix(),
'T': TGate().to_matrix(),
'0': np.array([[1, 0], [0, 0]], dtype=complex),
'1': np.array([[0, 0], [0, 1]], dtype=complex),
'+': np.array([[0.5, 0.5], [0.5, 0.5]], dtype=complex),
'-': np.array([[0.5, -0.5], [-0.5, 0.5]], dtype=complex),
'r': np.array([[0.5, -0.5j], [0.5j, 0.5]], dtype=complex),
'l': np.array([[0.5, 0.5j], [-0.5j, 0.5]], dtype=complex),
}
if re.match(r'^[IXYZHST01rl\-+]+$', label) is None:
raise QiskitError('Label contains invalid characters.')
# Initialize an identity matrix and apply each gate
num_qubits = len(label)
op = Operator(np.eye(2 ** num_qubits, dtype=complex))
for qubit, char in enumerate(reversed(label)):
if char != 'I':
op = op.compose(label_mats[char], qargs=[qubit])
return op
def is_unitary(self, atol=None, rtol=None):
"""Return True if operator is a unitary matrix."""
if atol is None:
atol = self.atol
if rtol is None:
rtol = self.rtol
return is_unitary_matrix(self._data, rtol=rtol, atol=atol)
def to_operator(self):
"""Convert operator to matrix operator class"""
return self
def to_instruction(self):
"""Convert to a UnitaryGate instruction."""
# pylint: disable=cyclic-import
from qiskit.extensions.unitary import UnitaryGate
return UnitaryGate(self.data)
def conjugate(self):
"""Return the conjugate of the operator."""
# Make a shallow copy and update array
ret = copy.copy(self)
ret._data = np.conj(self._data)
return ret
def transpose(self):
"""Return the transpose of the operator."""
# Make a shallow copy and update array
ret = copy.copy(self)
ret._data = np.transpose(self._data)
# Swap input and output dimensions
ret._set_dims(self._output_dims, self._input_dims)
return ret
def compose(self, other, qargs=None, front=False):
"""Return the composed operator.
Args:
other (Operator): an operator object.
qargs (list or None): a list of subsystem positions to apply
other on. If None apply on all
subsystems [default: None].
front (bool): If True compose using right operator multiplication,
instead of left multiplication [default: False].
Returns:
Operator: The operator self @ other.
Raise:
QiskitError: if operators have incompatible dimensions for
composition.
Additional Information:
Composition (``@``) is defined as `left` matrix multiplication for
matrix operators. That is that ``A @ B`` is equal to ``B * A``.
Setting ``front=True`` returns `right` matrix multiplication
``A * B`` and is equivalent to the :meth:`dot` method.
"""
if qargs is None:
qargs = getattr(other, 'qargs', None)
if not isinstance(other, Operator):
other = Operator(other)
# Validate dimensions are compatible and return the composed
# operator dimensions
input_dims, output_dims = self._get_compose_dims(
other, qargs, front)
# Full composition of operators
if qargs is None:
if front:
# Composition self * other
data = np.dot(self._data, other.data)
else:
# Composition other * self
data = np.dot(other.data, self._data)
return Operator(data, input_dims, output_dims)
# Compose with other on subsystem
if front:
num_indices = len(self._input_dims)
shift = len(self._output_dims)
right_mul = True
else:
num_indices = len(self._output_dims)
shift = 0
right_mul = False
# Reshape current matrix
# Note that we must reverse the subsystem dimension order as
# qubit 0 corresponds to the right-most position in the tensor
# product, which is the last tensor wire index.
tensor = np.reshape(self.data, self._shape)
mat = np.reshape(other.data, other._shape)
indices = [num_indices - 1 - qubit for qubit in qargs]
final_shape = [np.product(output_dims), np.product(input_dims)]
data = np.reshape(
Operator._einsum_matmul(tensor, mat, indices, shift, right_mul),
final_shape)
return Operator(data, input_dims, output_dims)
def dot(self, other, qargs=None):
"""Return the right multiplied operator self * other.
Args:
other (Operator): an operator object.
qargs (list or None): a list of subsystem positions to apply
other on. If None apply on all
subsystems [default: None].
Returns:
Operator: The operator self * other.
Raises:
QiskitError: if other cannot be converted to an Operator or has
incompatible dimensions.
"""
return super().dot(other, qargs=qargs)
def power(self, n):
"""Return the matrix power of the operator.
Args:
n (int): the power to raise the matrix to.
Returns:
BaseOperator: the n-times composed operator.
Raises:
QiskitError: if the input and output dimensions of the operator
are not equal, or the power is not a positive integer.
"""
if not isinstance(n, int):
raise QiskitError("Can only take integer powers of Operator.")
if self.input_dims() != self.output_dims():
raise QiskitError("Can only power with input_dims = output_dims.")
# Override base class power so we can implement more efficiently
# using Numpy.matrix_power
ret = copy.copy(self)
ret._data = np.linalg.matrix_power(self.data, n)
return ret
def tensor(self, other):
"""Return the tensor product operator self ⊗ other.
Args:
other (Operator): a operator subclass object.
Returns:
Operator: the tensor product operator self ⊗ other.
Raises:
QiskitError: if other cannot be converted to an operator.
"""
if not isinstance(other, Operator):
other = Operator(other)
input_dims = other.input_dims() + self.input_dims()
output_dims = other.output_dims() + self.output_dims()
data = np.kron(self._data, other._data)
return Operator(data, input_dims, output_dims)
def expand(self, other):
"""Return the tensor product operator other ⊗ self.
Args:
other (Operator): an operator object.
Returns:
Operator: the tensor product operator other ⊗ self.
Raises:
QiskitError: if other cannot be converted to an operator.
"""
if not isinstance(other, Operator):
other = Operator(other)
input_dims = self.input_dims() + other.input_dims()
output_dims = self.output_dims() + other.output_dims()
data = np.kron(other._data, self._data)
return Operator(data, input_dims, output_dims)
def _add(self, other, qargs=None):
"""Return the operator self + other.
If ``qargs`` are specified the other operator will be added
assuming it is identity on all other subsystems.
Args:
other (Operator): an operator object.
qargs (None or list): optional subsystems to add on
(Default: None)
Returns:
Operator: the operator self + other.
Raises:
QiskitError: if other is not an operator, or has incompatible
dimensions.
"""
# pylint: disable=import-outside-toplevel, cyclic-import
from qiskit.quantum_info.operators.scalar_op import ScalarOp
if qargs is None:
qargs = getattr(other, 'qargs', None)
if not isinstance(other, Operator):
other = Operator(other)
self._validate_add_dims(other, qargs)
other = ScalarOp._pad_with_identity(self, other, qargs)
ret = copy.copy(self)
ret._data = self.data + other.data
return ret
def _multiply(self, other):
"""Return the operator self * other.
Args:
other (complex): a complex number.
Returns:
Operator: the operator other * self.
Raises:
QiskitError: if other is not a valid complex number.
"""
if not isinstance(other, Number):
raise QiskitError("other is not a number")
ret = copy.copy(self)
ret._data = other * self._data
return ret
def equiv(self, other, rtol=None, atol=None):
"""Return True if operators are equivalent up to global phase.
Args:
other (Operator): an operator object.
rtol (float): relative tolerance value for comparison.
atol (float): absolute tolerance value for comparison.
Returns:
bool: True if operators are equivalent up to global phase.
"""
if not isinstance(other, Operator):
try:
other = Operator(other)
except QiskitError:
return False
if self.dim != other.dim:
return False
if atol is None:
atol = self.atol
if rtol is None:
rtol = self.rtol
return matrix_equal(self.data, other.data, ignore_phase=True,
rtol=rtol, atol=atol)
@property
def _shape(self):
"""Return the tensor shape of the matrix operator"""
return tuple(reversed(self.output_dims())) + tuple(
reversed(self.input_dims()))
@classmethod
def _einsum_matmul(cls, tensor, mat, indices, shift=0, right_mul=False):
"""Perform a contraction using Numpy.einsum
Args:
tensor (np.array): a vector or matrix reshaped to a rank-N tensor.
mat (np.array): a matrix reshaped to a rank-2M tensor.
indices (list): tensor indices to contract with mat.
shift (int): shift for indices of tensor to contract [Default: 0].
right_mul (bool): if True right multiply tensor by mat
(else left multiply) [Default: False].
Returns:
Numpy.ndarray: the matrix multiplied rank-N tensor.
Raises:
QiskitError: if mat is not an even rank tensor.
"""
rank = tensor.ndim
rank_mat = mat.ndim
if rank_mat % 2 != 0:
raise QiskitError(
"Contracted matrix must have an even number of indices.")
# Get einsum indices for tensor
indices_tensor = list(range(rank))
for j, index in enumerate(indices):
indices_tensor[index + shift] = rank + j
# Get einsum indices for mat
mat_contract = list(reversed(range(rank, rank + len(indices))))
mat_free = [index + shift for index in reversed(indices)]
if right_mul:
indices_mat = mat_contract + mat_free
else:
indices_mat = mat_free + mat_contract
return np.einsum(tensor, indices_tensor, mat, indices_mat)
@classmethod
def _init_instruction(cls, instruction):
"""Convert a QuantumCircuit or Instruction to an Operator."""
# Convert circuit to an instruction
if isinstance(instruction, QuantumCircuit):
instruction = instruction.to_instruction()
# Initialize an identity operator of the correct size of the circuit
op = Operator(np.eye(2 ** instruction.num_qubits))
op._append_instruction(instruction)
return op
@classmethod
def _instruction_to_matrix(cls, obj):
"""Return Operator for instruction if defined or None otherwise."""
if not isinstance(obj, Instruction):
raise QiskitError('Input is not an instruction.')
mat = None
if hasattr(obj, 'to_matrix'):
# If instruction is a gate first we see if it has a
# `to_matrix` definition and if so use that.
try:
mat = obj.to_matrix()
except QiskitError:
pass
return mat
def _append_instruction(self, obj, qargs=None):
"""Update the current Operator by apply an instruction."""
mat = self._instruction_to_matrix(obj)
if mat is not None:
# Perform the composition and inplace update the current state
# of the operator
op = self.compose(mat, qargs=qargs)
self._data = op.data
else:
# If the instruction doesn't have a matrix defined we use its
# circuit decomposition definition if it exists, otherwise we
# cannot compose this gate and raise an error.
if obj.definition is None:
raise QiskitError('Cannot apply Instruction: {}'.format(obj.name))
for instr, qregs, cregs in obj.definition:
if cregs:
raise QiskitError(
'Cannot apply instruction with classical registers: {}'.format(
instr.name))
# Get the integer position of the flat register
if qargs is None:
new_qargs = [tup.index for tup in qregs]
else:
new_qargs = [qargs[tup.index] for tup in qregs]
self._append_instruction(instr, qargs=new_qargs)