/
symplectic.py
2172 lines (1716 loc) · 76.4 KB
/
symplectic.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
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
Symplectic representation utility functions
"""
# ***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License. You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************
import numpy as _np
import copy as _copy
from pygsti.baseobjs.label import Label as _Label
from pygsti.baseobjs.smartcache import smart_cached
from pygsti.tools import matrixmod2 as _mtx
try:
from . import fastcalc as _fastcalc
except ImportError:
_fastcalc = None
def symplectic_form(n, convention='standard'):
"""
Creates the symplectic form for the number of qubits specified.
There are two variants, of the sympletic form over the finite field of the
integers modulo 2, used in pyGSTi. These corresponding to the 'standard' and
'directsum' conventions. In the case of 'standard', the symplectic form is the
2n x 2n matrix of ((0,1),(1,0)), where '1' and '0' are the identity and all-zeros
matrices of size n x n. The 'standard' symplectic form is probably the most
commonly used, and it is the definition used throughout most of the code,
including the Clifford compilers. In the case of 'directsum', the symplectic form
is the direct sum of n 2x2 bit-flip matrices. This is only used in pyGSTi for
sampling from the symplectic group.
Parameters
----------
n : int
The number of qubits the symplectic form should be constructed for. That
is, the function creates a 2n x 2n matrix that is a sympletic form
convention : str, optional
Can be either 'standard' or 'directsum', which correspond to two different
definitions for the symplectic form.
Returns
-------
numpy array
The specified symplectic form.
"""
nn = 2 * n
sym_form = _np.zeros((nn, nn), _np.int64)
assert(convention == 'standard' or convention == 'directsum')
if convention == 'standard':
sym_form[n:nn, 0:n] = _np.identity(n, _np.int64)
sym_form[0:n, n:nn] = _np.identity(n, _np.int64)
if convention == 'directsum':
# This current construction method is pretty stupid.
for j in range(0, n):
sym_form[2 * j, 2 * j + 1] = 1
sym_form[2 * j + 1, 2 * j] = 1
return sym_form
def change_symplectic_form_convention(s, outconvention='standard'):
"""
Maps the input symplectic matrix between the 'standard' and 'directsum' symplectic form conventions.
That is, if the input is a symplectic matrix with respect to the 'directsum'
convention and outconvention ='standard' the output of this function is the
equivalent symplectic matrix in the 'standard' symplectic form
convention. Similarily, if the input is a symplectic matrix with respect to the
'standard' convention and outconvention = 'directsum' the output of this function
is the equivalent symplectic matrix in the 'directsum' symplectic form
convention.
Parameters
----------
s : numpy.ndarray
The input symplectic matrix.
outconvention : str, optional
Can be either 'standard' or 'directsum', which correspond to two different
definitions for the symplectic form. This is the convention the input is
being converted to (and so the input should be a symplectic matrix in the
other convention).
Returns
-------
numpy array
The matrix `s` converted to `outconvention`.
"""
n = _np.shape(s)[0] // 2
if n == 1:
return _np.copy(s)
permutation_matrix = _np.zeros((2 * n, 2 * n), _np.int64)
for i in range(0, n):
permutation_matrix[2 * i, i] = 1
permutation_matrix[2 * i + 1, n + i] = 1
if outconvention == 'standard':
sout = _np.dot(_np.dot(permutation_matrix.T, s), permutation_matrix)
if outconvention == 'directsum':
sout = _np.dot(_np.dot(permutation_matrix, s), permutation_matrix.T)
return sout
def check_symplectic(m, convention='standard'):
"""
Checks whether a matrix is symplectic.
Parameters
----------
m : numpy array
The matrix to check.
convention : str, optional
Can be either 'standard' or 'directsum', Specifies the convention of
the symplectic form with respect to which the matrix should be
sympletic.
Returns
-------
bool
A bool specifying whether the matrix is symplectic
"""
n = _np.shape(m)[0] // 2
s_form = symplectic_form(n, convention=convention)
conj = _mtx.dot_mod2(_np.dot(m, s_form), _np.transpose(m))
return _np.array_equal(conj, s_form)
def inverse_symplectic(s):
"""
Returns the inverse of a symplectic matrix over the integers mod 2.
Parameters
----------
s : numpy array
The matrix to invert
Returns
-------
numpy array
The inverse of s, over the field of the integers mod 2.
"""
assert(check_symplectic(s)), "The input matrix is not symplectic!"
n = _np.shape(s)[0] // 2
s_form = symplectic_form(n)
s_inverse = _mtx.dot_mod2(_np.dot(s_form, _np.transpose(s)), s_form)
assert(check_symplectic(s_inverse)), "The inverse is not symplectic. Function has failed"
assert(_np.array_equal(_mtx.dot_mod2(s_inverse, s), _np.identity(2 * n, _np.int64))
), "The found matrix is not the inverse of the input. Function has failed"
return s_inverse
def inverse_clifford(s, p):
"""
Returns the inverse of a Clifford gate in the symplectic representation.
This uses the formualas derived in Hostens and De Moor PRA 71, 042315 (2005).
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
p : numpy array
The 'phase vector' over the integers mod 4 representing the Clifford
Returns
-------
sinverse : numpy array
The symplectic matrix representing the inverse of the input Clifford.
pinverse : numpy array
The 'phase vector' representing the inverse of the input Clifford.
"""
assert(check_valid_clifford(s, p)), \
"The input symplectic matrix - phase vector pair does not define a valid Clifford!"
sinverse = inverse_symplectic(s)
n = _np.shape(s)[0] // 2
# The formula used below for the inverse p vector comes from Hostens
# and De Moor PRA 71, 042315 (2005).
u = _np.zeros((2 * n, 2 * n), _np.int64)
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
vec1 = -1 * _np.dot(_np.transpose(sinverse), p)
inner = _np.dot(_np.dot(_np.transpose(sinverse), u), sinverse)
temp = 2 * _mtx.strictly_upper_triangle(inner) + _mtx.diagonal_as_matrix(inner)
temp = _mtx.diagonal_as_vec(_np.dot(_np.dot(_np.transpose(s), temp), s))
vec2 = -1 * _np.dot(_np.transpose(sinverse), temp)
vec3 = _mtx.diagonal_as_vec(inner)
pinverse = vec1 + vec2 + vec3
pinverse = pinverse % 4
assert(check_valid_clifford(sinverse, pinverse)), "The output does not define a valid Clifford. Function has failed"
s_check, p_check = compose_cliffords(s, p, sinverse, pinverse)
assert(_np.array_equal(s_check, _np.identity(2 * n, _np.int64))
), "The output is not the inverse of the input. Function has failed"
assert(_np.array_equal(p_check, _np.zeros(2 * n, _np.int64))
), "The output is not the inverse of the input. Function has failed"
s_check, p_check = compose_cliffords(sinverse, pinverse, s, p)
assert(_np.array_equal(s_check, _np.identity(2 * n, _np.int64))
), "The output is not the inverse of the input. Function has failed"
assert(_np.array_equal(p_check, _np.zeros(2 * n, _np.int64))
), "The output is not the inverse of the input. Function has failed"
return sinverse, pinverse
def check_valid_clifford(s, p):
"""
Checks if a symplectic matrix - phase vector pair (s,p) is the symplectic representation of a Clifford.
This uses the formualas derived in Hostens and De Moor PRA 71, 042315 (2005).
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
p : numpy array
The 'phase vector' over the integers mod 4 representing the Clifford
Returns
-------
bool
True if (s,p) is the symplectic representation of some Clifford.
"""
# Checks if the matrix s is symplectic, which is the only constraint on s.
is_symplectic_matrix = check_symplectic(s)
# Check whether the phase vector is valid. This currently does *not* check
# that p is a vector over [0,1,2,3]. Perhaps it should. The constraint
# that we check is satisfied comes from Hostens and De Moor PRA 71, 042315 (2005).
n = _np.shape(s)[0] // 2
u = _np.zeros((2 * n, 2 * n), _np.int64)
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
vec = p + _mtx.diagonal_as_vec(_np.dot(_np.dot(_np.transpose(s), u), s))
vec = vec % 2
is_valid_phase_vector = _np.array_equal(vec, _np.zeros(len(p), _np.int64))
assert(is_valid_phase_vector)
return (is_symplectic_matrix and is_valid_phase_vector)
def construct_valid_phase_vector(s, pseed):
"""
Constructs a phase vector that, when paired with the provided symplectic matrix, defines a Clifford gate.
If the seed phase vector, when paired with `s`, represents some Clifford this
seed is returned. Otherwise 1 mod 4 is added to the required elements of the
`pseed` in order to make it at valid phase vector (which is one of many possible
phase vectors that, together with s, define a valid Clifford).
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
pseed : numpy array
The seed 'phase vector' over the integers mod 4.
Returns
-------
numpy array
Some p such that (s,p) is the symplectic representation of some Clifford.
"""
pout = pseed.copy()
n = _np.shape(s)[0] // 2
assert(check_symplectic(s)), "The input matrix is not symplectic!"
u = _np.zeros((2 * n, 2 * n), _np.int64)
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
# Each element of this vector should be 0 (mod 2) if this is a valid phase vector.
# This comes from the formulas in Hostens and De Moor PRA 71, 042315 (2005).
vec = pout + _mtx.diagonal_as_vec(_np.dot(_np.dot(_np.transpose(s), u), s))
vec = vec % 2
# Adds 1 mod 4 to all the elements of the vector where the required constraint is
# not satisfied. This is then always a valid phase vector.
pout[vec != 0] += 1
pout = pout % 4
assert(check_valid_clifford(s, pout)), "The output does not define a valid Clifford. Function has failed"
return pout
def find_postmultipled_pauli(s, p_implemented, p_target, qubit_labels=None):
"""
Finds the Pauli layer that should be appended to a circuit to implement a given Clifford.
If some circuit implements the clifford described by the symplectic matrix s and
the vector p_implemented, this function returns the Pauli layer that should be
appended to this circuit to implement the clifford described by s and the vector
p_target.
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
implemented by the circuit
p_implemented : numpy array
The 'phase vector' over the integers mod 4 representing the Clifford
implemented by the circuit
p_target : numpy array
The 'phase vector' over the integers mod 4 that, together with `s` represents
the Clifford that you want to implement. Together with `s`, this vector must
define a valid Clifford.
qubit_labels : list, optional
A list of qubit labels, that are strings or ints. The length of this list should
be equal to the number of qubits the Clifford acts on. The ith element of the
list is the label corresponding to the qubit at the ith index of `s` and the two
phase vectors. If None, defaults to the integers from 0 to number of qubits - 1.
Returns
-------
list
A list that defines a Pauli layer, with the ith element containig one of the
4 tuples (P,qubit_labels[i]) with P = 'I', 'Z', 'Y' and 'Z'
"""
n = _np.shape(s)[0] // 2
s_form = symplectic_form(n)
vec = _mtx.dot_mod2(s, _np.dot(s_form, (p_target - p_implemented) // 2))
if qubit_labels is None:
qubit_labels = list(range(n))
pauli_layer = []
for q in range(0, n):
if vec[q] == 0 and vec[q + n] == 0:
pauli_layer.append(('I', qubit_labels[q]))
elif vec[q] == 0 and vec[q + n] == 1:
pauli_layer.append(('Z', qubit_labels[q]))
elif vec[q] == 1 and vec[q + n] == 0:
pauli_layer.append(('X', qubit_labels[q]))
elif vec[q] == 1 and vec[q + n] == 1:
pauli_layer.append(('Y', qubit_labels[q]))
return pauli_layer
def find_premultipled_pauli(s, p_implemented, p_target, qubit_labels=None):
"""
Finds the Pauli layer that should be prepended to a circuit to implement a given Clifford.
If some circuit implements the clifford described by the symplectic matrix s and
the vector p_implemented, this function returns the Pauli layer that should be
prefixed to this circuit to implement the clifford described by s and the vector
p_target.
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
implemented by the circuit
p_implemented : numpy array
The 'phase vector' over the integers mod 4 representing the Clifford
implemented by the circuit
p_target : numpy array
The 'phase vector' over the integers mod 4 that, together with `s` represents
the Clifford that you want to implement. Together with `s`, this vector must
define a valid Clifford.
qubit_labels : list, optional
A list of qubit labels, that are strings or ints. The length of this list should
be equal to the number of qubits the Clifford acts on. The ith element of the
list is the label corresponding to the qubit at the ith index of `s` and the two
phase vectors. If None, defaults to the integers from 0 to number of qubits - 1.
Returns
-------
list
A list that defines a Pauli layer, with the ith element containig one of the
4 tuples ('I',i), ('X',i), ('Y',i), ('Z',i).
"""
n = _np.shape(s)[0] // 2
s_form = symplectic_form(n)
vec = _mtx.dot_mod2(s_form, (p_target - p_implemented) // 2)
if qubit_labels is None:
qubit_labels = list(range(n))
pauli_layer = []
for q in range(n):
if vec[q] == 0 and vec[q + n] == 0:
pauli_layer.append(('I', qubit_labels[q]))
elif vec[q] == 0 and vec[q + n] == 1:
pauli_layer.append(('Z', qubit_labels[q]))
elif vec[q] == 1 and vec[q + n] == 0:
pauli_layer.append(('X', qubit_labels[q]))
elif vec[q] == 1 and vec[q + n] == 1:
pauli_layer.append(('Y', qubit_labels[q]))
return pauli_layer
def find_pauli_layer(pvec, qubit_labels, pauli_labels=None):
"""
TODO: docstring
pauli_labels defaults to ['I', 'X', 'Y', 'Z'].
"""
if pauli_labels is None:
pauli_labels = ['I', 'X', 'Y', 'Z']
paulis_as_int_list = find_pauli_number(pvec)
return [(pauli_labels[p], q) for p, q in zip(paulis_as_int_list, qubit_labels)]
def find_pauli_number(pvec):
"""
TODO: docstring
"""
n = len(pvec) // 2
v = (pvec[0:n] // 2) + 2 * (pvec[n:] // 2)
return [[0, 3, 1, 2][i] for i in v] # [0,0]=I, [2,0]=Z, [0,2]=X, and [2,2]=Y.
def compose_cliffords(s1, p1, s2, p2, do_checks=True):
"""
Multiplies two cliffords in the symplectic representation.
The output corresponds to the symplectic representation of C2 times C1 (i.e., C1
acts first) where s1 (s2) and p1 (p2) are the symplectic matrix and phase vector,
respectively, for Clifford C1 (C2). This uses the formualas derived in Hostens
and De Moor PRA 71, 042315 (2005).
Parameters
----------
s1 : numpy array
The symplectic matrix over the integers mod 2 representing the first Clifford
p1 : numpy array
The 'phase vector' over the integers mod 4 representing the first Clifford
s2 : numpy array
The symplectic matrix over the integers mod 2 representing the second Clifford
p2 : numpy array
The 'phase vector' over the integers mod 4 representing the second Clifford
do_checks : bool
If True (default), check inputs and output are valid cliffords.
If False, these checks are skipped (for speed)
Returns
-------
s : numpy array
The symplectic matrix over the integers mod 2 representing the composite Clifford
p : numpy array
The 'phase vector' over the integers mod 4 representing the compsite Clifford
"""
assert(_np.shape(s1) == _np.shape(s2)), "Input must be Cliffords acting on the same number of qubits!"
if do_checks:
assert(check_valid_clifford(s1, p1)), "The first matrix-vector pair is not a valid Clifford!"
assert(check_valid_clifford(s2, p2)), "The second matrix-vector pair is not a valid Clifford!"
n = _np.shape(s1)[0] // 2
# Below we calculate the s and p for the composite Clifford using the formulas from
# Hostens and De Moor PRA 71, 042315 (2005).
s = _mtx.dot_mod2(s2, s1)
u = _np.zeros((2 * n, 2 * n), _np.int64)
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
vec1 = _np.dot(_np.transpose(s1), p2)
inner = _np.dot(_np.dot(_np.transpose(s2), u), s2)
matrix = 2 * _mtx.strictly_upper_triangle(inner) + _mtx.diagonal_as_matrix(inner)
vec2 = _mtx.diagonal_as_vec(_np.dot(_np.dot(_np.transpose(s1), matrix), s1))
vec3 = _np.dot(_np.transpose(s1), _mtx.diagonal_as_vec(inner))
p = p1 + vec1 + vec2 - vec3
p = p % 4
if do_checks:
assert(check_valid_clifford(s, p)), "The output is not a valid Clifford! Function has failed."
return s, p
def symplectic_kronecker(sp_factors):
"""
Takes a kronecker product of symplectic representations.
Construct a single `(s,p)` symplectic (or stabilizer) representation that
corresponds to the tensor (kronecker) product of the objects represented
by each `(s,p)` element of `sp_factors`.
This is performed by inserting each factor's `s` and `p` elements into the
appropriate places of the final (large) `s` and `p` arrays. This operation
works for combining Clifford operations AND also stabilizer states.
Parameters
----------
sp_factors : iterable
A list of `(s,p)` symplectic (or stabilizer) representation factors.
Returns
-------
s : numpy.ndarray
An array of shape (2n,2n) where n is the *total* number of qubits (the
sum of the number of qubits in each `sp_factors` element).
p : numpy.ndarray
A 1D array of length 2n.
"""
nlist = [len(p) // 2 for s, p in sp_factors] # number of qubits per factor
n = sum(nlist) # total number of qubits
sout = _np.zeros((2 * n, 2 * n), _np.int64)
pout = _np.zeros(2 * n, _np.int64)
k = 0 # current qubit index
for (s, p), nq in zip(sp_factors, nlist):
assert(s.shape == (2 * nq, 2 * nq))
sout[k:k + nq, k:k + nq] = s[0:nq, 0:nq]
sout[k:k + nq, n + k:n + k + nq] = s[0:nq, nq:2 * nq]
sout[n + k:n + k + nq, k:k + nq] = s[nq:2 * nq, 0:nq]
sout[n + k:n + k + nq, n + k:n + k + nq] = s[nq:2 * nq, nq:2 * nq]
pout[k:k + nq] = p[0:nq]
pout[n + k:n + k + nq] = p[nq:2 * nq]
k += nq
return sout, pout
def prep_stabilizer_state(nqubits, zvals=None):
"""
Contruct the `(s,p)` stabilizer representation for a computational basis state given by `zvals`.
Parameters
----------
nqubits : int
Number of qubits
zvals : iterable, optional
An iterable over anything that can be cast as True/False
to indicate the 0/1 value of each qubit in the Z basis.
If None, the all-zeros state is created. If None, then
all zeros is assumed.
Returns
-------
s,p : numpy.ndarray
The stabilizer "matrix" and phase vector corresponding to the desired
state. `s` has shape (2n,2n) (it includes antistabilizers) and `p`
has shape 2n, where n equals `nqubits`.
"""
n = nqubits
s = _np.fliplr(_np.identity(2 * n, _np.int64)) # flip b/c stab cols are *first*
p = _np.zeros(2 * n, _np.int64)
if zvals:
for i, z in enumerate(zvals):
p[i] = p[i + n] = 2 if bool(z) else 0 # EGN TODO: check this is right -- (how to update the destabilizers?)
return s, p
def apply_clifford_to_stabilizer_state(s, p, state_s, state_p):
"""
Applies a clifford in the symplectic representation to a stabilizer state in the standard stabilizer representation.
The output corresponds to the stabilizer representation of the output state.
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
p : numpy array
The 'phase vector' over the integers mod 4 representing the Clifford
state_s : numpy array
The matrix over the integers mod 2 representing the stabilizer state
state_p : numpy array
The 'phase vector' over the integers mod 4 representing the stabilizer state
Returns
-------
out_s : numpy array
The symplectic matrix over the integers mod 2 representing the output state
out_p : numpy array
The 'phase vector' over the integers mod 4 representing the output state
"""
two_n = _np.shape(s)[0]; n = two_n // 2
assert(_np.shape(state_s) == (two_n, two_n)), "Clifford and state must be for the same number of qubits!"
assert(_np.shape(state_p) == (two_n,)), "Invalid stabilizer state representation"
assert(check_valid_clifford(s, p)), "The `s`,`p` matrix-vector pair is not a valid Clifford!"
#EGN TODO: check valid stabilizer state?
# Below we calculate the s and p for the output state using the formulas from
# Hostens and De Moor PRA 71, 042315 (2005).
out_s = _mtx.dot_mod2(s, state_s)
u = _np.zeros((2 * n, 2 * n), _np.int64)
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
inner = _np.dot(_np.dot(_np.transpose(s), u), s)
vec1 = _np.dot(_np.transpose(state_s), p - _mtx.diagonal_as_vec(inner))
matrix = 2 * _mtx.strictly_upper_triangle(inner) + _mtx.diagonal_as_matrix(inner)
vec2 = _mtx.diagonal_as_vec(_np.dot(_np.dot(_np.transpose(state_s), matrix), state_s))
out_p = state_p + vec1 + vec2
out_p = out_p % 4
##More explicitly operates on stabilizer and antistabilizer separately, but same as above
#out_p = _np.zeros(2*n, _np.int64)
#for slc in (slice(0,n),slice(n,2*n)):
# ss = state_s[:,slc]
# inner = _np.dot(_np.dot(_np.transpose(s),u),s)
# vec1 = _np.dot(_np.transpose(ss),p - _mtx.diagonal_as_vec(inner))
# matrix = 2*_mtx.strictly_upper_triangle(inner)+_mtx.diagonal_as_matrix(inner)
# vec2 = _mtx.diagonal_as_vec(_np.dot(_np.dot(_np.transpose(ss),matrix),ss))
# out_p[slc] = state_p[slc] + vec1 + vec2
# out_p[slc] = out_p[slc] % 4
#EGN TODO: check for valid stabilizer state
return out_s, out_p
def pauli_z_measurement(state_s, state_p, qubit_index):
"""
Computes the probabilities of 0/1 (+/-) outcomes from measuring a Pauli operator on a stabilizer state.
Parameters
----------
state_s : numpy array
The matrix over the integers mod 2 representing the stabilizer state
state_p : numpy array
The 'phase vector' over the integers mod 4 representing the stabilizer state
qubit_index : int
The index of the qubit being measured
Returns
-------
p0, p1 : float
Probabilities of 0 (+ eigenvalue) and 1 (- eigenvalue) outcomes.
state_s_0, state_s_1 : numpy array
Matrix over the integers mod 2 representing the output stabilizer
states.
state_p_0, state_p_1 : numpy array
Phase vectors over the integers mod 4 representing the output
stabilizer states.
"""
two_n = len(state_p); n = two_n // 2
assert(_np.shape(state_s) == (two_n, two_n)), "Inconsistent stabilizier representation!"
#This algorithm follows that of PRA 70, 052328 (2004),
# except note:
# 0) columns & rows are reversed
# 1) that states carry arount full "mod 4"
# phase vectors instead of the minmal "mod 2" ones
# comprised of the higher-order-bit of each "mod 4" vector.
# 2) that the stabilizer is held in the *first* n
# columns of state_s rather than the last n.
a = qubit_index
#print("DB: pauli Z_%d meas on state:" % a) #DEBUG
#print(state_s); print(state_p) #DEBUG
# Let A be the Pauli op being measured
#Step1: determine if all stabilizer elements commute with Z_a
# which amounts to checking whether there are any 1-bits in the a-th
# row of state_s for the first n columns (any stabilizers that have X's
# or Y's for qubit a, as 00=I, 10=X, 11=Y, 01=Z)
for col in range(n):
if state_s[a, col] == 1:
p = col
#print("Column %d anticommutes with Z_%d" % (p,a)) #DEBUG
# p is first stabilizer that anticommutes w/Z_a. Outcome is random,
# and we just need to update the state (if requested).
s_out = state_s.copy(); p_out = state_p.copy()
for i in range(two_n):
if i != p and state_s[a, i] == 1:
#print("COLSUM ",i,p) #DEBUG
colsum(i, p, s_out, p_out, n)
#print(s_out); print(p_out) #DEBUG
#print("-----") #DEBUG
s_out[:, p + n] = s_out[:, p]; p_out[p + n] = p_out[p] # set p-th col -> (p+n)-th col
s_out[:, p] = 0; s_out[a + n, p] = 1 # p-th col = I*...Z_a*...I stabilizer
icount = sum([3 if (s_out[i, p] == s_out[i + n, p] == 1) else 0 for i in range(n)]) # 11 = -iY convention
p_out0 = p_out.copy(); p_out0[p] = (4 - (icount % 4)) % 4 # so overall phase is 0
p_out1 = p_out.copy(); p_out1[p] = (p_out0[p] + 2) % 4 # so overall phase is -1
return 0.5, 0.5, s_out, s_out, p_out0, p_out1 # Note: s is same for 0 and 1 outcomes
#print("Nothing anticommutes!") #DEBUG
# no break ==> all commute, so outcome is deterministic, so no
# state update; just determine whether Z_a or -Z_a is in the stabilizer,
# which we do using the "anti-stabilizer" cleverness of PRA 70, 052328
acc_s = _np.zeros(two_n, _np.int64); acc_p = _np.zeros(1, _np.int64)
for i in range(n, two_n): # loop over anti-stabilizer
if state_s[a, i] == 1: # for elements that anti-commute w/Z_a
colsum_acc(acc_s, acc_p, i - n, state_s, state_p, n) # act w/corresponding *stabilizer* el
# now the high bit of acc_p holds the outcome
icount = acc_p[0] + sum([3 if (acc_s[i] == acc_s[i + n] == 1) else 0 for i in range(n)]) # 11 = -iY convention
icount = icount % 4
if icount == 0: # outcome is always zero
#print("Always 0!") #DEBUG
return (1.0, 0.0, state_s, state_s, state_p, state_p)
else: # outcome is always 1
#print("Always 1!") #DEBUG
assert(icount == 2) # should never get 1 or 3 (low bit should always be 0)
return (0.0, 1.0, state_s, state_s, state_p, state_p)
def colsum(i, j, s, p, n):
"""
A helper routine used for manipulating stabilizer state representations.
Updates the `i`-th stabilizer generator (column of `s` and element of `p`)
with the group-action product of the `j`-th and the `i`-th generators, i.e.
generator[i] -> generator[j] + generator[i]
Parameters
----------
i : int
Destination generator index.
j : int
Sournce generator index.
s : numpy array
The matrix over the integers mod 2 representing the stabilizer state
p : numpy array
The 'phase vector' over the integers mod 4 representing the stabilizer state
n : int
The number of qubits. `s` must be shape (2n,2n) and `p` must be
length 2n.
Returns
-------
None
"""
#OLD: according to Aaronson convention (not what we use)
##Note: need to divide p-vals by 2 to access high-bit
#test = 2*(p[i]//2 + p[j]//2) + sum(
# [colsum_g(s[k,j], s[k+n,j], s[k,i], s[k+n,i]) for k in range(n)] )
#test = test % 4
#assert(test in (0,2)) # test should never be congruent to 1 or 3 (mod 4)
#p[i] = 0 if (test == 0) else 2 # ( = 10 = 1 in high bit)
u = _np.zeros((2 * n, 2 * n), _np.int64) # CACHE!
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
p[i] += p[j] + 2 * float(_np.dot(s[:, i].T, _np.dot(u, s[:, j])))
for k in range(n):
s[k, i] = s[k, j] ^ s[k, i]
s[k + n, i] = s[k + n, j] ^ s[k + n, i]
#EGN TODO: use _np.bitwise_xor or logical_xor here? -- keep it obvious (&slow) for now...
return
def colsum_acc(acc_s, acc_p, j, s, p, n):
"""
A helper routine used for manipulating stabilizer state representations.
Similar to :func:`colsum` except a separate "accumulator" column is
used instead of the `i`-th column of `s` and element of `p`. I.e., this
performs:
acc[0] -> generator[j] + acc[0]
Parameters
----------
acc_s : numpy array
The matrix over the integers mod 2 representing the "accumulator" stabilizer state
acc_p : numpy array
The 'phase vector' over the integers mod 4 representing the "accumulator" stabilizer state
j : int
Index of the stabilizer generator being accumulated (see above).
s : numpy array
The matrix over the integers mod 2 representing the stabilizer state
p : numpy array
The 'phase vector' over the integers mod 4 representing the stabilizer state
n : int
The number of qubits. `s` must be shape (2n,2n) and `p` must be
length 2n.
Returns
-------
None
"""
##Note: need to divide p-vals by 2 to access high-bit
#test = 2*(acc_p[0]//2 + p[j]//2) + sum(
# [colsum_g(s[k,j], s[k+n,j], acc_s[k], acc_s[k+n]) for k in range(n)] )
#test = test % 4
#assert(test in (0,2)) # test should never be congruent to 1 or 3 (mod 4)
#acc_p[0] = 0 if (test == 0) else 2 # ( = 10 = 1 in high bit)
u = _np.zeros((2 * n, 2 * n), _np.int64) # CACHE!
u[n:2 * n, 0:n] = _np.identity(n, _np.int64)
acc_p[0] += p[j] + 2 * float(_np.dot(acc_s.T, _np.dot(u, s[:, j])))
for k in range(n):
acc_s[k] = s[k, j] ^ acc_s[k]
acc_s[k + n] = s[k + n, j] ^ acc_s[k + n]
#EGN TODO: use _np.bitwise_xor or logical_xor here? -- keep it obvious (&slow) for now...
return
def stabilizer_measurement_prob(state_sp_tuple, moutcomes, qubit_filter=None,
return_state=False):
"""
Compute the probability of a given outcome when measuring some or all of the qubits in a stabilizer state.
Returns this probability, optionally along with the updated (post-measurement) stabilizer state.
Parameters
----------
state_sp_tuple : tuple
A `(s,p)` tuple giving the stabilizer state to measure.
moutcomes : array-like
The z-values identifying which measurement outcome (a computational
basis state) to compute the probability for.
qubit_filter : iterable, optional
If not None, a list of qubit indices which are measured.
`len(qubit_filter)` should always equal `len(moutcomes)`. If None, then
assume *all* qubits are measured (`len(moutcomes)` == num_qubits).
return_state : bool, optional
Whether the post-measurement (w/outcome `moutcomes`) state is also
returned.
Returns
-------
p : float
The probability of the given measurement outcome.
state_s,state_p : numpy.ndarray
Only returned when `return_state=True`. The post-measurement stabilizer
state representation (an updated version of `state_sp_tuple`).
"""
state_s, state_p = state_sp_tuple # should be a StabilizerState.to_dense() "object"
p = 1
if qubit_filter is None: # len(moutcomes) == nQubits
qubit_filter = range(len(moutcomes))
for i, outcm in zip(qubit_filter, moutcomes):
p0, p1, ss0, ss1, sp0, sp1 = pauli_z_measurement(state_s, state_p, i)
# could cache these results in a FUTURE _stabilizer_measurement_probs function?
if outcm == 0:
p *= p0; state_s, state_p = ss0, sp0
else:
p *= p1; state_s, state_p = ss1, sp1
return (p, state_s, state_p) if return_state else p
def embed_clifford(s, p, qubit_inds, n):
"""
Embeds the `(s,p)` Clifford symplectic representation into a larger symplectic representation.
The action of `(s,p)` takes place on the qubit indices specified by `qubit_inds`.
Parameters
----------
s : numpy array
The symplectic matrix over the integers mod 2 representing the Clifford
p : numpy array
The 'phase vector' over the integers mod 4 representing the Clifford
qubit_inds : list
A list or array of integers specifying which qubits `s` and `p` act on.
n : int
The total number of qubits
Returns
-------
s : numpy array
The symplectic matrix over the integers mod 2 representing the embedded Clifford
p : numpy array
The 'phase vector' over the integers mod 4 representing the embedded Clifford
"""
ne = len(qubit_inds) # nQubits for embedded_op
s_out = _np.identity(2 * n, _np.int64)
p_out = _np.zeros(2 * n, _np.int64)
for i, di in enumerate(qubit_inds): # di = "destination index"
p_out[di] = p[i]
p_out[di + n] = p[i + ne]
for j, dj in enumerate(qubit_inds):
s_out[di, dj] = s[i, j]
s_out[di + n, dj + n] = s[i + ne, j + ne]
s_out[di, dj + n] = s[i, j + ne]
s_out[di + n, dj] = s[i + ne, j]
return s_out, p_out
def compute_internal_gate_symplectic_representations(gllist=None):
"""
Creates a dictionary of the symplectic representations of 'standard' Clifford gates.
Returns a dictionary containing the symplectic matrices and phase vectors that represent
the specified 'standard' Clifford gates, or the representations of *all* the standard gates
if no list of operation labels is supplied. These 'standard' Clifford gates are those gates that
are already known to the code gates (e.g., the label 'CNOT' has a specfic meaning in the
code), and are recorded as unitaries in "internalgates.py".
Parameters
----------
gllist : list, optional
If not None, a list of strings corresponding to operation labels for any of the standard
gates that have fixed meaning for the code (e.g., 'CNOT' corresponds to
the CNOT gate with the first qubit the target). For example, this list could be
gllist = ['CNOT','H','P','I','X'].
Returns
-------
srep_dict : dict
dictionary of `(smatrix,svector)` tuples, where `smatrix` and `svector`
are numpy arrays containing the symplectic matrix and phase vector
representing the operation label given by the key.
"""
# Full dictionaries, containing the symplectic representations of *all* gates
# that are hard-coded, and which have a specific meaning to the code.
complete_s_dict = {}
complete_p_dict = {}
# The Pauli gates
complete_s_dict['I'] = _np.array([[1, 0], [0, 1]], _np.int64)
complete_s_dict['X'] = _np.array([[1, 0], [0, 1]], _np.int64)
complete_s_dict['Y'] = _np.array([[1, 0], [0, 1]], _np.int64)
complete_s_dict['Z'] = _np.array([[1, 0], [0, 1]], _np.int64)
complete_p_dict['I'] = _np.array([0, 0], _np.int64)
complete_p_dict['X'] = _np.array([0, 2], _np.int64)
complete_p_dict['Y'] = _np.array([2, 2], _np.int64)
complete_p_dict['Z'] = _np.array([2, 0], _np.int64)
# Five single qubit gates that each represent one of five classes of Cliffords
# that equivalent up to Pauli gates and are not equivalent to idle (that class
# is covered by any one of the Pauli gates above).
complete_s_dict['H'] = _np.array([[0, 1], [1, 0]], _np.int64)
complete_s_dict['P'] = _np.array([[1, 0], [1, 1]], _np.int64)
complete_s_dict['PH'] = _np.array([[0, 1], [1, 1]], _np.int64)
complete_s_dict['HP'] = _np.array([[1, 1], [1, 0]], _np.int64)
complete_s_dict['HPH'] = _np.array([[1, 1], [0, 1]], _np.int64)
complete_p_dict['H'] = _np.array([0, 0], _np.int64)
complete_p_dict['P'] = _np.array([1, 0], _np.int64)
complete_p_dict['PH'] = _np.array([0, 1], _np.int64)
complete_p_dict['HP'] = _np.array([3, 0], _np.int64)
complete_p_dict['HPH'] = _np.array([0, 3], _np.int64)
# The full 1-qubit Cliffor group, using the same labelling as in extras.rb.group
complete_s_dict['C0'] = _np.array([[1, 0], [0, 1]], _np.int64)
complete_p_dict['C0'] = _np.array([0, 0], _np.int64)
complete_s_dict['C1'] = _np.array([[1, 1], [1, 0]], _np.int64)
complete_p_dict['C1'] = _np.array([1, 0], _np.int64)
complete_s_dict['C2'] = _np.array([[0, 1], [1, 1]], _np.int64)
complete_p_dict['C2'] = _np.array([0, 1], _np.int64)