Skip to content

Commit

Permalink
alternative equivalence checking algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
@sebastiaanbrand
sebastiaanbrand committed Sep 30, 2024
1 parent 10e37d3 commit 6ec786c
Show file tree
Hide file tree
Showing 3 changed files with 142 additions and 9 deletions.
93 changes: 90 additions & 3 deletions examples/circuit_equivalence.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ static size_t max_tablesize = 1LL<<25;
static size_t min_cachesize = 1LL<<16;
static size_t max_cachesize = 1LL<<18;
static size_t wgt_tab_size = 1LL<<23;
static char eqcheck_alg[20] = "alternating";
static double tolerance = 1e-14;
static double threshold = 1e-6;
static int wgt_table_type = COMP_HASHMAP;
Expand All @@ -32,7 +33,8 @@ static struct argp_option options[] =
{"workers", 'w', "<workers>", 0, "Number of workers/threads (default=1)", 0},
{"norm-strat", 's', "<low|max|min|l2>", 0, "Edge weight normalization strategy (default=max)", 0},
{"tol", 't', "<tolerance>", 0, "Tolerance for deciding edge weights equal (default=1e-14)", 0},
{"threshold", 1000, "<threshold>", 0, "Threshold for deciding when two matrices are close enough to be equivalent (default=1e-6)", 0},
{"algorithm", 1000, "<pauli|alternating>", 0, "Which equivalence checking algorithm to use.", 0},
{"threshold", 1001, "<threshold>", 0, "Threshold for deciding when two matrices are close enough to be equivalent (default=1e-6)", 0},
{"count-nodes", 'c', 0, 0, "Track maximum number of nodes", 0},
{0, 0, 0, 0, 0, 0}
};
Expand All @@ -57,6 +59,11 @@ parse_opt(int key, char *arg, struct argp_state *state)
count_nodes = true;
break;
case 1000:
if (strcmp(arg, "pauli") != 0 && strcmp(arg, "alternating") != 0)
argp_usage(state);
strncpy(eqcheck_alg, arg, 20);
break;
case 1001:
threshold = atof(arg);
break;
case ARGP_KEY_ARG:
Expand Down Expand Up @@ -97,6 +104,7 @@ void print_stats() {
// print stats in JSON format
printf("{\n");
printf(" \"statistics\": {\n");
printf(" \"algorithm\": \"%s\",\n", eqcheck_alg);
printf(" \"circuit_U\": \"%s\",\n", circuit_U->name);
printf(" \"circuit_V\": \"%s\",\n", circuit_V->name);
printf(" \"counterexample\": \"%s\",\n", stats.counterexample);
Expand Down Expand Up @@ -186,6 +194,7 @@ QMDD get_gate_matrix(quantum_op_t* gate, BDDVAR nqubits, bool dag) {
}
}


QMDD compute_UPUdag(quantum_circuit_t *circuit, gate_id_t P, BDDVAR k) {
BDDVAR nqubits = circuit->qreg_size;
QMDD circ_matrix = qmdd_create_single_qubit_gate(nqubits, k, P);
Expand Down Expand Up @@ -218,7 +227,13 @@ QMDD compute_UPUdag(quantum_circuit_t *circuit, gate_id_t P, BDDVAR k) {
return circ_matrix;
}

void circuit_equivalence_check(quantum_circuit_t *U, quantum_circuit_t *V) {

/**
* Equivalence checking based on Theorem 1 of
* "Fast equivalence checking of quantum circuits of Clifford gates",
* D. Thanos et al., ATVA (2023)
*/
void pauli_echeck(quantum_circuit_t *U, quantum_circuit_t *V) {
if (U->qreg_size != V->qreg_size) {
snprintf(stats.counterexample, sizeof(stats.counterexample),
"Circuits have different number of qubits");
Expand Down Expand Up @@ -266,6 +281,73 @@ QMDD compute_UPUdag(quantum_circuit_t *circuit, gate_id_t P, BDDVAR k) {
}
}


/**
* Equivalence checking using the "alternating" algorith from
* "Advanced Equivalence Checking for Quantum Circuits",
* L. Burgholzer, R. Wille, (TCAD), 2021
*/
void alternating_eqcheck(quantum_circuit_t *U, quantum_circuit_t *V)
{
// check number of qubits
if (U->qreg_size != V->qreg_size) {
snprintf(stats.counterexample, sizeof(stats.counterexample),
"Circuits have different number of qubits");
return;
}
BDDVAR nqubits = U->qreg_size;

// get circuits as arrays + set U as longest circuit
int ngates_u, ngates_v;
quantum_op_t **ops_u = circuit_as_array(U, &ngates_u);
quantum_op_t **ops_v = circuit_as_array(V, &ngates_v);
if (ngates_u < ngates_v) {
quantum_op_t **tmp1 = ops_v; ops_v = ops_u; ops_u = tmp1;
int tmp2 = ngates_v; ngates_v = ngates_u; ngates_u = tmp2;
}

// compute V^dagger * U from the inside out,
// in steps of 1 for V, and steps of |U|/|V| for U
double step_u = (double)ngates_u / (double)ngates_v;
int i = 0;
double j = 0;
QMDD prod, vi_dag, uj;
prod = qmdd_create_all_identity_matrix(nqubits);
while (i < ngates_v) {
// compute V[i]^dagger * prod * U[j-step_u] * ... * U[j]
vi_dag = get_gate_matrix(ops_v[i], nqubits, true);
prod = evbdd_matmat_mult(vi_dag, prod, nqubits);
for (int _j = ceilf(j); (_j < j + step_u) && (_j < ngates_u) ; _j += 1) {
uj = get_gate_matrix(ops_u[_j], nqubits, false);
prod = evbdd_matmat_mult(prod, uj, nqubits);
}
i += 1;
j += step_u;
}
assert (i == ngates_v);
assert (ceilf(j) == ngates_u);

// check if results is equal to identity matrix
QMDD identity = qmdd_create_all_identity_matrix(nqubits);
if (prod == identity) {
stats.fidelity = 1;
stats.equivalent = true;
}
else {
// not exactly identity, check fidelity based on given threshold
double norm = pow(2.0, nqubits*2);
double fid = qmdd_fidelity(prod, identity, nqubits*2) / norm;
stats.fidelity = fid;
if (fabs(fid - 1.0) < threshold)
stats.equivalent = true;
else
stats.equivalent = false;
}

free(ops_u);
free(ops_v);
}

int main(int argc, char *argv[])
{
argp_parse(&argp, argc, argv, 0, 0, 0);
Expand All @@ -279,7 +361,12 @@ int main(int argc, char *argv[])
clock_t cpu_t1 = clock();
double wall_t1 = wctime();

circuit_equivalence_check(circuit_U, circuit_V);
if (strcmp(eqcheck_alg, "pauli") == 0)
pauli_echeck(circuit_U, circuit_V);
else if (strcmp(eqcheck_alg, "alternating") == 0)
alternating_eqcheck(circuit_U, circuit_V);
else
fprintf(stderr, "Invalid arg for algorithm (should be caught by argp)\n");

clock_t cpu_t2 = clock();
double wall_t2 = wctime();
Expand Down
41 changes: 39 additions & 2 deletions qasm/qsylvan_qasm_parser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -715,8 +715,45 @@ void optimize_qubit_order(quantum_circuit_t *circuit, bool allow_swaps)
}
}

quantum_op_t** circuit_as_array(quantum_circuit_t *circuit, int *length)
{
// loop over circuit to get lenght
*length = 0;
quantum_op_t *op = circuit->operations;
while (op != NULL) {
switch (op->type) {
case op_gate:
case op_measurement:
*length += 1;
break;
default:
break;
}
op = op->next;
}

// loop over circuit and put into array
quantum_op_t **ops_array = (quantum_op_t**)malloc(sizeof(quantum_op_t*) * (*length));
int i = 0;
op = circuit->operations;
while (op != NULL) {
switch (op->type) {
case op_gate:
case op_measurement:
ops_array[i] = op;
i++;
break;
default:
break;
}
op = op->next;
}

return ops_array;
}


void print_quantum_circuit(quantum_circuit_t* circuit)
void print_quantum_circuit(quantum_circuit_t *circuit)
{
printf("qreg size: %d\n", circuit->qreg_size);
printf("creg size: %d\n", circuit->creg_size);
Expand All @@ -728,7 +765,7 @@ void print_quantum_circuit(quantum_circuit_t* circuit)
}
}

void fprint_creg(FILE *stream, quantum_circuit_t* circuit)
void fprint_creg(FILE *stream, quantum_circuit_t *circuit)
{
// print in big-endian
for (int i = circuit->creg_size-1; i >= 0 ; i--) {
Expand Down
17 changes: 13 additions & 4 deletions qasm/qsylvan_qasm_parser.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,16 +74,25 @@ void optimize_qubit_order(quantum_circuit_t *circuit, bool allow_swaps);


/**
* Free all quantum elements found in quantum_op_s including 'first'. What is 'first'?
* Free all quantum elements found in quantum_op_s including 'first'.
*/
void free_quantum_circuit(quantum_circuit_t* circuit);
void free_quantum_circuit(quantum_circuit_t *circuit);

/**
* Return the circuit as array of quantum_op_t's instead of linked list.
* 'length' is used to return the length of the array (number of operations).
*
* The returned array is malloced and should be freed with free().
* The quantum_circuit_t itself should be freed separately with free_quantum_circuit().
*/
quantum_op_t** circuit_as_array(quantum_circuit_t *circuit, int *length);

/**
* Print operations of quantum circuit, quantum elements or control registers.
*/
void print_quantum_op(quantum_op_t *op);
void print_quantum_circuit(quantum_circuit_t* circuit);
void fprint_creg(FILE *stream, quantum_circuit_t* circuit);
void print_quantum_circuit(quantum_circuit_t *circuit);
void fprint_creg(FILE *stream, quantum_circuit_t *circuit);


#ifdef __cplusplus
Expand Down

0 comments on commit 6ec786c

Please sign in to comment.