Quantum Computing Simulation

Outline

• Simulation: what and why

• Classical and Reversible Computing

• Classical Bits

• Quantum Bits

• Measurements

• Entanglement

• Algorithms - addition, Grover, QFT

• Simulation = software implementation of a mathematical model

• Why? Study algorithms, experiment with parameters and noise

• Advantages: see all state amplitudes, probabilities, entanglement

• Disadvantage: need O(2ⁿ) space and time for n qubits

  Simulation limit: n ≈ 30 on single node, n ≈ 50 distributed on 65000 nodes

• Pure states (state vector) vs. mixed states (density matrix) vs. unentangled (separable)

Classical Computing Separate input and output bits Generally not reversible Distributed in space		Reversible Computing Input bits become the output bits Distributed in time
Ion Trap Quantum Device

One Classical Bit

4 possible functions

~   x  f0 f1 f2 f3     reversible: f1 = x,  f2 = NOT x
~   0   0  0  1  1
~   1   0  1  0  1

Two Classical Bits

16 possible functions

~  x y  f0 f1 f2 f3 f4 f5 f6 f7 f8 f9 fa fb fc fd fe ff
~  0 0   0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
~  0 1   0  0  0  0  1  1  1  1  0  0  0  0  1  1  1  1
~  1 0   0  0  1  1  0  0  1  1  0  0  1  1  0  0  1  1
~  1 1   0  1  0  1  0  1  0  1  0  1  0  1  0  1  0  1

~ (x,y) -> (x,f(x,y)) reversible: f5 = y,  fa = NOT y,  f6 = x XOR y,  f9 = NOT (x XOR y)

~ f6: CNOT(x,y) -> (x, x XOR y)

Simulating One Classical Bit

void X( int b[]) { b[0] = !b[0]; } // NOT

int main( void) { int b[1] = { 0 }; // constraint: b[i] = 0 or 1
...

Simulating Two (or more) Classical Bits

void X( int b[], int k) { b[k] = !b[k]; } // NOT on bit k

void CX( int b[], int c, int k) // controlled-NOT on bit k controlled by bit c, c != k
{
~  b[k] ^= b[c]; // if( b[c]) b[k] = !b[k];
}

int main( void) { int b[2] = { 0, 0 }; 
...

Simulating One Quantum Bit

typedef double complex Complex;

// swap a[i] and a[j]
void swap( Complex a[], int i, int j) { Complex t = a[i]; a[i] = a[j]; a[j] = t; }

void X( Complex a[]) { swap( a, 0, 1); } // NOT

void Z( Complex a[]) { a[1] = -a[1]; } // phase flip

void S( Complex a[]) { a[1] = CMPLX(0,1)*a[1]; } // sqrt(Z)

void H( Complex a[]) // Hadamard
{ 
~  double s2 = 1/sqrt(2);

~  Complex t0 = a[0], t1 = a[1]; a[0] = s2*(t0+t1); a[1] = s2*(t0-t1);
}
...
int main( void) { Complex a[2] = { 3.0/5.0, CMPLX(0,4.0/5.0) }; // constraint: sum |a[i]|2 = 1
...

Single Qubit General Unitary Operation

// OpenQASM3 U with global phase gamma
//
void U( Complex a[], double theta, double phi, double lambda, double gamma)
{
~  double t = theta/2; Complex t0 = a[0], t1 = a[1], p = exp(CMPLX(0,gamma)),

~  c = p*cos(t), s = p*sin(t), u00 = c, u01 = -exp(CMPLX(0,lambda))*s,

~  u10 = exp(CMPLX(0,phi))*s, u11 = exp(CMPLX(0,phi+lambda))*c;

~  a[0] = u00*t0+u01*t1; a[1] = u10*t0+u11*t1;
}

Simulating Two Quantum Bits

//    q1 q0
// a0  0 0
// a1  0 1
// a2  1 0
// a3  1 1
//
void X( Complex a[], int k) // NOT on qubit k
{
~  if( k == 0) { swap( a, 0, 1); swap( a, 2, 3); }

~  else        { swap( a, 0, 2); swap( a, 1, 3); }
}
void CX( Complex a[], int c, int k) // controlled-NOT on qubit k controlled by qubit c
{
~  if( c == 1) swap( a, 2, 3); else swap( a, 1, 3);
}
void Z( Complex a[], int k) // phase-flip qubit k
{
~  if( k == 0) { a[1] = -a[1]; a[3] = -a[3]; }

~  else        { a[2] = -a[2]; a[3] = -a[3]; }
}
...
int main( void) { Complex a[4] = { 1, 0, 0, 0 };
...

Simulating Two Quantum Bits - Hadamard

//    q1 q0
// a0  0 0
// a1  0 1
// a2  1 0
// a3  1 1
//
void h( Complex a[], int i, int j) // half Hadamard
{
~  Complex t0 = a[i], t1 = a[j]; a[i] = s2*(t0+t1); a[j] = s2*(t0-t1);
}
void H( Complex a[], int k) // Hadamard on qubit k
{
~  if( k == 0) { h( a, 0, 1); h( a, 2, 3); }

~  else        { h( a, 0, 2); h( a, 1, 3); }
}

Simulating n Quantum Bits - Hadamard

// Hadamard on qubit k
//
// n qubits, N = 2n amplitudes
//
~ unsigned long i = 0, j, k1 = 1LU << k; // k1 = 0...010...0, k'th bit = 1

~ while( i < N) { // i will have k'th bit = 0, j will have k'th bit = 1

~   j = i | k1; Complex t0 = a[i], t1 = a[j];

~   a[i] = s2*(t0+t1); a[j] = s2*(t0-t1);

~   ++i; i += (i & k1); // skip i values with k'th bit = 1

~ }

Simulating Measurements

One Qubit

// complex magnitude squared
//
double abs2( Complex c) { double x = creal(c), y = cimag(c); return x*x + y*y; }

//  0.0  ... r ...       1.0
//   |<------->|<-------->|
//       p0         p1
// 
void M( Complex a[]) { // measure

~  double r = rand() / (RAND_MAX + 1.0); // 0.0 <= r < 1.0

~  double p = abs2( a[0]);

~  if( r < p) { printf( "M = 0\\\\n"); a[0] = 1; a[1] = 0; }

~  else       { printf( "M = 1\\\\n"); a[0] = 0; a[1] = 1; }
}

Simulating Measurements

Two Qubits - measure both

//
//  0.0       ... r ...                       1.0
//   |<------->|<-------->|<------->|<-------->|
//       p00        p01       p10       p11
//
int Mk( Complex a[]) { // measure overall state

~  double r = rand() / (RAND_MAX + 1.0); // 0.0 <= r < 1.0

~  double p = 0; // cumulative probability

~  int k = 3; for( int i = 0; i < 4; ++i) { p += abs2( a[i]); if( r < p) { k = i; break; } }

~  printf( "k = %i\\\\n", k); return k;
}

void M10( Complex a[]) { // measure qubits 1 and 0

~  int k = Mk( a); a[0]=a[1]=a[2]=a[3]=0; a[k] = 1;
}

Simulating Measurements

Two Qubits - measure one

void M1( Complex a[]) { // measure qubit 1

~  int k = Mk( a); if( k == 0 || k == 1) { // qubit 1 is 0

~    printf( "M1 = 0\\\\n"); a[2] = a[3] = 0; // normalize remaining states

~    double s = sqrt( abs2(a[0]) + abs2(a[1])); a[0] /= s; a[1] /= s; }

~  else { // qubit 1 is 1

~    printf( "M1 = 1\\\\n"); a[0] = a[1] = 0;

~    double s = sqrt( abs2(a[2]) + abs2(a[3])); a[2] /= s; a[3] /= s; }
}

Measuring Entanglement - Theory

 

• trace distance = 0 for no entanglement, 1 for max entanglement = 1-fidelity

   

• fidelity = (|<A,B>|/(||A||*||B||))²

  where A (B) is the reduced state with qubit k equal to 0 (1)

   

• If qubit k is not entangled with any of the other qubits,

  then the reduced state A should be the same as the reduced state B,

  i.e. the reduced state should be independent of qubit k.

   

• <A,B>/(||A||*||B||)) is the cosine of the angle between vectors A and B.

   

• A and B represent the same state if the angle is 0 or 180 degrees,

  which corresponds to the cosine having a magnitude squared value of 1.

Measuring Entanglement - Code

Real tangle( unsigned int k) const {

~ Complex prod = 0; Real A = 0, B = 0, r = 0;

~ unsigned long i = 0, j, k1 = 1LU << k; // k1 = 0...010...0, k'th bit = 1

~ while( i < N) { // i will have k'th bit = 0, j will have k'th bit = 1

~   j = i | k1; prod += Conj(a[i])*a[j];

~   A += abs2(a[i]); B += abs2(a[j]);

~   ++i; i += (i & k1); // skip i values with k'th bit = 1
~ }

~ Real t = 2*eps*eps*N; // threshold for zero norm squared: 2*sum{|e+i*e|**2} = 2*(2*e*e*(N/2))

~ if( A > t && B > t) r = 1 - abs2(prod)/(A*B);

~ return r;
}

Entanglement Example - Super Dense Coding

• Alice and Bob start with two qubits in state |00> with random global phase.

• After performing H1, then CX10, Bob takes qubit #0 and goes away.

  E is the qubit entanglement measure.

  ~    q: (0.0849528,0.996385) (0,0) (0,0) (0,0)
  ~    E: 0 0
  ~   H1: (0.0600707,0.704551) (0,0) (0.0600707,0.704551) (0,0)
  ~    E: 0 0
  ~ CX10: (0.0600707,0.704551) (0,0) (0,0) (0.0600707,0.704551)
  ~    E: 1 1
  

• Later, to communicate the value 3, Alice performs X1, then Z1, then sends qubit #1 to Bob.

  ~   X1: (0,0) (0.0600707,0.704551) (0.0600707,0.704551) (0,0)
  ~    E: 1 1
  ~   Z1: (0,0) (0.0600707,0.704551) (-0.0600707,-0.704551) (0,0)
  ~    E: 1 1
  

• Bob then performs CX10, then H1, then measures.

  ~ CX10: (0,0) (0.0600707,0.704551) (0,0) (-0.0600707,-0.704551)
  ~    E: 0 0
  ~   H1: (0,0) (0,0) (0,0) (0.0849528,0.996385)
  ~    E: 0 0
  ~  P(3) = 1
  

Super Dense Coding with Hadamard Phase Shift Error (p/p₀ = 0.8)

Alice + Bob

~    q: (0.285804,0.958288) (0,0) (0,0) (0,0)
~    E: 0 0
~   H1: (0.202094,0.677612) (0,0) (-0.0171906,0.706898) (0,0)
~    E: 0 0
~ CX10: (0.202094,0.677612) (0,0) (0,0) (-0.0171906,0.706898)
~    E: 1 1

Alice

~   X1: (0,0) (-0.0171906,0.706898) (0.202094,0.677612) (0,0)
~    E: 1 1
~   Z1: (0,0) (-0.0171906,0.706898) (-0.202094,-0.677612) (0,0)
~    E: 1 1

Bob

~ CX10: (0,0) (-0.0171906,0.706898) (0,0) (-0.202094,-0.677612)
~    E: 0 0
~   H1: (0,0) (-0.296127,0.0883184) (0,0) (-0.0231213,0.950775)
~    E: 0 0
~  P(3) = 0.904508

Addition

low-level implementation

Addition

high-level simulation, (x,y)->(x,(x+y)%M)

    I = state index

    x = I >> m;       // top m bits, M = 2m

    y = I & mask;     // bottom m bits, mask = M-1 = 0111..1 (m 1's)

    j = (x + y) % M;  // classical addition

    J = (x << m) | j; // concatenate the bits

    tj = a[J];        // temporarily save original a[J]

    a[J] = a[I];      // update a[J]

Grover's Algorithm

theory

Grover's Algorithm

low-level implementation

Grover's Algorithm

high-level simulation

// flip sign of state m
//
q.a[m] = -q.a[m]; 

// inversion about the average
//
Complex avg = 0; for( unsigned long i = 0; i < q.N; ++i) avg += q.a[i];

avg *= 2.0/q.N; for( unsigned long i = 0; i < q.N; ++i) q.a[i] = avg - q.a[i];

Grover's Algorithm

Entanglement

Quantum Fourier Transform

low-level implementation

Quantum Fourier Transform

high-level simulation

~ // quantum FFT
~ //
~ void qfft( int inv = 0) {

~   fft<Complex,Real>( a, N, inv); // classical FFT, without scaling and bit reversal

~  // scale values and bit-reverse the indices
~  //
~   Real s = 1/std::sqrt((Real)N);

~   for( unsigned long i = 0, j = 0; i < N; ++i) {

~     j = bitrev( i, n);

~     if( i < j) { swap( i, j); a[i] *= s; a[j] *= s; }

~     else if( i == j) a[i] *= s;
~   }
~ }

Future Work

 

• Simulated quantum state is the same,

  simulated at a low-level or simulated at a high-level.

  Do both and compare.

   

• Is simulated quantum state the same as physical implementation?

  Probably not, depending on noise.

   

• How to realistically simulate noise?

   

• How to design quantum algorithms resistant to noise?