Sequential Quantum Gate Decomposer  v1.9.3
Powerful decomposition of general unitarias into one- and two-qubit gates gates
Random_Unitary.cpp
Go to the documentation of this file.
1 /*
2 Created on Fri Jun 26 14:13:26 2020
3 Copyright 2020 Peter Rakyta, Ph.D.
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 
17 You should have received a copy of the GNU General Public License
18 along with this program. If not, see http://www.gnu.org/licenses/.
19 
20 @author: Peter Rakyta, Ph.D.
21 */
27 #include "Random_Unitary.h"
28 
29 
36 Matrix
37 few_CNOT_unitary( int qbit_num, int cnot_num) {
38 
39  // the current number of CNOT gates
40  int cnot_num_curr = 0;
41 
42  // the size of the matrix
43  int matrix_size = Power_of_2(qbit_num);
44 
45  // The unitary describing each qubits in their initial state
46  Matrix mtx = create_identity( matrix_size );
47 
48  // constructing the unitary
49  while (true) {
50  int cnot_or_u3 = rand() % 5 + 1;
51 
52  CNOT* cnot_op = NULL;
53  U3* u3_op = NULL;
54 
55  Matrix gate_matrix;
56 
57  if (cnot_or_u3 <= 4) {
58  // creating random parameters for the U3 operation
59  Matrix_real parameters(1, 3);
60  parameters[0] = double(rand())/RAND_MAX*4*M_PI;
61  parameters[1] = double(rand())/RAND_MAX*2*M_PI;
62  parameters[2] = double(rand())/RAND_MAX*2*M_PI;
63 
64 
65  // randomly choose the target qbit
66  int target_qbit = rand() % qbit_num;
67 
68  // creating the U3 gate
69  u3_op = new U3(qbit_num, target_qbit);
70 
71  // get the matrix of the operation
72 
73  gate_matrix = u3_op->get_matrix(parameters);
74  }
75  else if ( cnot_or_u3 == 5 ) {
76  // randomly choose the target qbit
77  int target_qbit = rand() % qbit_num;
78 
79  // randomly choose the control qbit
80  int control_qbit = rand() % qbit_num;
81 
82  if (target_qbit == control_qbit) {
83  gate_matrix = create_identity( matrix_size );
84  }
85  else {
86 
87  // creating the CNOT gate
88  cnot_op = new CNOT(qbit_num, control_qbit, target_qbit);
89 
90  // get the matrix of the operation
91  gate_matrix = cnot_op->get_matrix();
92 
93  cnot_num_curr = cnot_num_curr + 1;
94  }
95  }
96  else {
97  gate_matrix = create_identity( matrix_size );
98  }
99 
100 
101  // get the current unitary
102  Matrix mtx_tmp = dot(gate_matrix, mtx);
103  mtx = mtx_tmp;
104 
105  delete u3_op;
106  u3_op = NULL;
107 
108  delete cnot_op;
109  cnot_op = NULL;
110 
111  // exit the loop if the maximal number of CNOT gates reached
112  if (cnot_num_curr >= cnot_num) {
113  return mtx;
114  }
115 
116  }
117 
118 }
119 
120 
127 
128  if (dim_in < 2) {
129  throw "wrong dimension";
130  }
131 
132  // number of qubits
133  dim = dim_in;
134 
135 }
136 
137 
142 Matrix
144 
145  // create array of random parameters to construct random unitary
146  double* vartheta = (double*) qgd_calloc( int(dim*(dim-1)/2),sizeof(double), 64);
147  double* varphi = (double*) qgd_calloc( int(dim*(dim-1)/2),sizeof(double), 64);
148  double* varkappa = (double*) qgd_calloc( (dim-1),sizeof(double), 64);
149 
150  // initialize random seed:
151  srand (time(NULL));
152 
153  for (int idx=0; idx<dim*(dim-1)/2; idx++) {
154  vartheta[idx] = (2*double(rand())/double(RAND_MAX)-1)*2*M_PI;
155  }
156 
157 
158  for (int idx=0; idx<dim*(dim-1)/2; idx++) {
159  varphi[idx] = (2*double(rand())/double(RAND_MAX)-1)*2*M_PI;
160  }
161 
162 
163  for (int idx=0; idx<(dim-1); idx++) {
164  varkappa[idx] = (2*double(rand())/double(RAND_MAX)-1)*2*M_PI;
165  }
166 
167  Matrix Umtx = Construct_Unitary_Matrix( vartheta, varphi, varkappa );
168 
169  qgd_free(vartheta);
170  qgd_free(varphi);
171  qgd_free(varkappa);
172  vartheta = NULL;
173  varphi = NULL;
174  varkappa = NULL;
175 
176  return Umtx;
177 
178 }
179 
180 
188 Matrix
189 Random_Unitary::Construct_Unitary_Matrix( double* vartheta, double* varphi, double* varkappa ) {
190 
191 
192  Matrix ret = create_identity(dim);
193 
194  for (int varalpha=1; varalpha<dim; varalpha++) { // = 2:obj.dim
195  for (int varbeta = 0; varbeta<varalpha; varbeta++) {// 1:varalpha-1
196 
197  double theta_loc = vartheta[ convert_indexes(varalpha, varbeta) ];
198  double phi_loc = varphi[ convert_indexes(varalpha, varbeta) ];
199 
200 
201  // Eq (26)
202  QGD_Complex16 a;
203  a.real = cos( theta_loc )*cos(phi_loc);
204  a.imag = cos( theta_loc )*sin(phi_loc);
205 
206  // Eq (28) and (26)
207  double varepsilon = varkappa[varalpha-1]*kronecker( varalpha-1, varbeta);
208  QGD_Complex16 b;
209  b.real = sin( theta_loc )*cos(varepsilon);
210  b.imag = sin( theta_loc )*sin(varepsilon);
211 
212  a.real = -a.real;
213  b.imag = -b.imag;
214  Matrix Omega_loc = Omega( varalpha, varbeta, a, b );
215  Matrix ret_tmp = dot( ret, Omega_loc); // ret * Omega_loc
216 
217  ret = ret_tmp;
218  }
219  }
220 
221 
222  QGD_Complex16 gamma_loc;
223  gamma_loc.real = gamma();
224  gamma_loc.imag = 0;
225 
226  for ( int idx=0; idx<dim*dim; idx++) {
227  ret[idx] = mult(ret[idx], gamma_loc);
228  }
229 
230  return ret;
231 
232 
233 
234 }
235 
236 
243 int
244 Random_Unitary::convert_indexes( int varalpha, int varbeta ) {
245  int ret = varbeta + (varalpha-1)*(varalpha-2)/2;
246  return ret;
247 }
248 
249 
256  return Construct_Unitary_Matrix( parameters, parameters+int(dim*(dim-1)/2), parameters+int(dim*(dim-1)));
257 }
258 
259 
268 Matrix
269 Random_Unitary::Omega(int varalpha, int varbeta, QGD_Complex16 x, QGD_Complex16 y ) {
270 
271  Matrix ret = I_alpha_beta( varalpha, varbeta );
272 
273 
274  Matrix Mloc;
275 
276  if (varalpha + varbeta != (3 + kronecker( dim, 2 )) ) {
277  Mloc = M( varalpha, varbeta, x, y );
278 
279  }
280  else {
281  y.imag = -y.imag;
282  Mloc = M( varalpha, varbeta, x, mult(gamma(), y) );
283  }
284 
285 
286  //#pragma omp parallel for
287  for ( int idx=0; idx<dim*dim; idx++ ) {
288  ret[idx].real = ret[idx].real + Mloc[idx].real;
289  ret[idx].imag = ret[idx].imag + Mloc[idx].imag;
290  }
291 
292  return ret;
293 
294 
295 }
296 
297 
306 Matrix
307 Random_Unitary::M( int varalpha, int varbeta, QGD_Complex16 s, QGD_Complex16 t ) {
308 
309  Matrix Qloc = Q( s, t);
310 
311  Matrix ret1 = E_alpha_beta( varbeta, varbeta );
312  Matrix ret2 = E_alpha_beta( varbeta, varalpha );
313  Matrix ret3 = E_alpha_beta( varalpha, varbeta );
314  Matrix ret4 = E_alpha_beta( varalpha, varalpha );
315 
316 
317  mult(Qloc[0], ret1);
318  mult(Qloc[1], ret2);
319  mult(Qloc[2], ret3);
320  mult(Qloc[3], ret4);
321 
322  Matrix ret = Matrix(dim, dim);
323 
324  for ( int idx=0; idx<dim*dim; idx++ ) {
325  ret[idx].real = ret1[idx].real + ret2[idx].real + ret3[idx].real + ret4[idx].real;
326  ret[idx].imag = ret1[idx].imag + ret2[idx].imag + ret3[idx].imag + ret4[idx].imag;
327  }
328 
329  return ret;
330 
331 }
332 
333 
340 Matrix
342 
343  Matrix ret = Matrix(2, 2);
344  ret[0] = u2;
345  ret[1] = u1;
346  ret[2].real = -u1.real;
347  ret[2].imag = u1.imag;
348  ret[3].real = u2.real;
349  ret[3].imag = -u2.imag;
350 
351  return ret;
352 
353 }
354 
355 
362 Matrix
363 Random_Unitary::E_alpha_beta( int varalpha, int varbeta ) {
364 
365  Matrix ret = Matrix(dim, dim);
366  memset( ret.get_data(), 0, dim*dim*sizeof(QGD_Complex16));
367  ret[varalpha*dim+varbeta].real = 1;
368 
369  return ret;
370 
371 }
372 
379 Matrix
380 Random_Unitary::I_alpha_beta( int varalpha, int varbeta ) {
381 
382 
383  Matrix ret = create_identity(dim);
384 
385  ret[varalpha*dim+varalpha].real = 0;
386  ret[varbeta*dim+varbeta].real = 0;
387 
388  return ret;
389 
390 }
391 
392 
397 double
399 
400  double ret = pow(-1, 0.25*(2*dim-1+pow(-1,dim)));//(-1)^(0.25*(2*dim-1+(-1)^dim));
401  return ret;
402 
403 }
404 
411 double
412 Random_Unitary::kronecker( int a, int b ) {
413 
414  if (a == b) {
415  return 1;
416  }
417  else {
418  return 0;
419  }
420 
421 }
422 
423 
424 
425 
426 
427 
Matrix dot(Matrix &A, Matrix &B)
Call to calculate the product of two complex matrices by calling method zgemm3m from the CBLAS librar...
Definition: dot.cpp:38
A class representing a U3 gate.
Definition: U3.h:36
Matrix few_CNOT_unitary(int qbit_num, int cnot_num)
Call to create a random unitary constructed by CNOT operation between randomly chosen qubits and by r...
void * qgd_calloc(int element_num, int size, int alignment)
custom defined memory allocation function.
Definition: common.cpp:74
Matrix get_matrix(Matrix_real &parameters)
Call to retrieve the gate matrix.
Definition: U3.cpp:102
void qgd_free(void *ptr)
custom defined memory release function.
Definition: common.cpp:102
QGD_Complex16 mult(QGD_Complex16 &a, QGD_Complex16 &b)
Call to calculate the product of two complex scalars.
Definition: common.cpp:259
int convert_indexes(int varalpha, int varbeta)
Calculates an index from paramaters varalpha and varbeta.
double kronecker(int a, int b)
Kronecker delta.
scalar * get_data() const
Call to get the pointer to the stored data.
int dim
The number of rows in the created unitary.
double gamma()
Implements Eq (11) of arXiv:1303:5904v1.
matrix_size
[load Umtx]
Definition: example.py:58
Matrix E_alpha_beta(int varalpha, int varbeta)
Implements matrix I below Eq (7) of arXiv:1303:5904v1.
Umtx
The unitary to be decomposed.
Definition: example.py:53
Matrix M(int varalpha, int varbeta, QGD_Complex16 s, QGD_Complex16 t)
Implements Eq (8) of arXiv:1303:5904v1.
Structure type representing complex numbers in the SQUANDER package.
Definition: QGDTypes.h:38
A class representing a CNOT operation.
Definition: CNOT.h:35
int Power_of_2(int n)
Calculates the n-th power of 2.
Definition: common.cpp:117
Class to store data of complex arrays and its properties.
Definition: matrix.h:38
Random_Unitary(int dim_in)
Constructor of the class.
Matrix Q(QGD_Complex16 u1, QGD_Complex16 u2)
Implements Eq (9) of arXiv:1303:5904v1.
Matrix I_alpha_beta(int varalpha, int varbeta)
Implements matrix I below Eq (7) of arXiv:1303:5904v1.
Matrix create_identity(int matrix_size)
Call to create an identity matrix.
Definition: common.cpp:164
double real
the real part of a complex number
Definition: QGDTypes.h:40
Matrix Omega(int varalpha, int varbeta, QGD_Complex16 x, QGD_Complex16 y)
Eq (6) of arXiv:1303:5904v1.
Matrix get_matrix()
Call to retrieve the operation matrix.
Definition: CNOT.cpp:118
Matrix Construct_Unitary_Matrix()
Call to create a random unitary.
Class to store data of complex arrays and its properties.
Definition: matrix_real.h:39
double imag
the imaginary part of a complex number
Definition: QGDTypes.h:42