Quantum Computing Simulator

Introduction

There is a ‘new’ repository on GitHub¹. It’s not that new, I think I started it in September 2022, but I didn’t add any description on it here. So, here it is now: A quantum computation project. It’s a quantum computing simulator together with many algorithms as examples and tests. It started first as a very simple simulator, one that simply followed the math, with matrices, tensor products and so on, without any enhancement, but I got carried away, I added more and more example algorithms and it got quite slow… besides, I couldn’t implement some things (that needed oracles and so on, needing quite a bit of ancilla qubits) so I changed the implementation switching from the old algorithm that was O(N^2) to one that is O(N) (where N is the number of basis states, which increases exponentially with the number of qubits). That allowed me to implement more algorithms and to re-implement some already implemented, but fully, not providing the oracles/functions just as big matrices as I did before. This optimization is a big deal, at the moment I’ve done it, it went from well over 5 minutes of runtime to something I think lower than 20 seconds. Now with the added algorithms, the optimized one runs on my computer in about 20 seconds, the one non-optimized… well, I don’t think I have patience to wait for it to finish execution.

Algorithms

Here are the example algorithms, as a list:

Grover
Deutsch-Jozsa
Simon
Quantum Fourier Transform
Phase estimation
Shor
Bernstein–Vazirani
Quantum counting
Quantum teleportation
Entanglement swapping
Superdense coding
Quantum cryptography: BB84 protocol
CHSH inequality violation

Quantum error correction:

3-qubit error correcting a qubit-flip
3-qubit error correcting a sign-flip
Shor Code

Quantum adders:

Quantum half-adder for 1-qubit
Quantum full-adder for 1-qubit
Full adder for two N-qubit numbers
Draper adder
Draper adder with carry

Simulation of quantum simulation:

Evolution in time with a Hamiltonian given as a sum of Pauli products with real coefficients
Evolution of a 1D Gaussian packet in time, solving the 1D time-dependent Schrodinger equation using a Trotter decomposition and quantum Fourier transform

Paradoxes (although some of the above might be considered paradoxes as well):

Quantum eraser
General Elitzur-Vaidman Bomb tester/interaction free measurement/counterfactual computation
Hardy’s paradox

Quantum games:

Coin flipping
Quantum pseudo-telepathy: The magic square game

Distributed quantum computing:

Distributed CNOT (derived from a class that allows distribution of any 2-qubits controlled gate)
Quantum CNOT gate teleportation

Quantum Machine Learning

2D Q-Means Clustering

Optimization

QAOA (quantum approximate optimization algorithm) on Ising Model
VQE (variational quantum eigensolver)

Explaining them would probably take an entire book, I don’t have enough patience for that, so I can recommend a book instead: Quantum Computation and Quantum Information²

Theory

I won’t present the theory here, except give some links to some of it. The basics are quite simple and the already mentioned book is very accessible.
Here are a couple of papers that I used for some of the examples: Introduction to Coding Quantum Algorithms: A Tutorial Series Using Qiskit³, Fundamentals In Quantum Algorithms: A Tutorial Series Using Qiskit Continued⁴. This one has many algoritms, it’s good also for ideas to try out some other ones, I didn’t implement them all: Quantum Algorithm Implementations for Beginners⁵. An easy one allowing you to start implementing your own simulator and Grover and Shor algorithms: Undergraduate computational physics projects on quantum computing⁶.

The code

First, as many other open source projects I have on GitHub, the project uses Eigen⁷.

The code for the simulator is quite simple, although the whole thing is more complex due of the various algorithms implemented. Because there is quite a bit of it, I’ll try to point out the most important parts.
The most important classes are in QubitRegister.h and QuantumGate.h. Those headers alone (well, together with the ones they include, most notably Eigen and QubitRegisterCalculator.h) are enough for quantum computing simulation.
Here is the minimal code to run a simulation:

#include <iostream>

#include "QubitRegister.h"

int main()
{
    QC::QubitRegister reg(3);
    QC::Gates::PauliXGate x;
    reg.ApplyGate(x, 0);
    unsigned int m = reg.MeasureAll();

    std::cout << "Result should be 1, it's: " << m << std::endl;

    return 0;
}

#include <iostream>

#include "QubitRegister.h"

int main()

{

QC::QubitRegister reg(3);

QC::Gates::PauliXGate x;

reg.ApplyGate(x, 0);

unsigned int m = reg.MeasureAll();

std::cout << "Result should be 1, it's: " << m << std::endl;

return 0;

}

In the QuantumGate.h, there is a QuantumGate class which is no big deal, it’s just a base class as a starting point. QuantumGateWithOp is just adding some implementation and has the Eigen matrix for the gate operator as a member. Then SingleQubitGate, TwoQubitsGate, and ThreeQubitsGate follow. Those are to be used as base classes for specific gates implementations. Here is the class for Hadamard:

template<class MatrixClass = Eigen::MatrixXcd> class HadamardGate : public SingleQubitGate<MatrixClass>
{
public:
    using BaseClass = SingleQubitGate<MatrixClass>;
    using OpClass = typename BaseClass::BaseClass;

    HadamardGate()
    {
        // 1./sqrt(2) * (X + Z)
        static const double norm = 1. / sqrt(2.);
        OpClass::operatorMat(0, 0) = norm;
        OpClass::operatorMat(0, 1) = norm;
        OpClass::operatorMat(1, 0) = norm;
        OpClass::operatorMat(1, 1) = -norm;
    }
};

template<class MatrixClass = Eigen::MatrixXcd> class HadamardGate : public SingleQubitGate<MatrixClass>

{

public:

using BaseClass = SingleQubitGate<MatrixClass>;

using OpClass = typename BaseClass::BaseClass;

HadamardGate()

{

// 1./sqrt(2) * (X + Z)

static const double norm = 1. / sqrt(2.);

OpClass::operatorMat(0, 0) = norm;

OpClass::operatorMat(0, 1) = norm;

OpClass::operatorMat(1, 0) = norm;

OpClass::operatorMat(1, 1) = -norm;

}

};

As most projects presenting quantum computing simulation use the method with matrices, I preserved the old code in there (commenting out #define OPTIMIZED_TENSOR_PRODUCT 1 in QubitRegister.h falls back to it, but I strongly advise against doing that). I think it’s easier to understand by approaching the subject that way. Here is the function that calculates the tensor product for a single qubit gate (the matrix that is going to multiply the ‘register’ vector, to apply the operator on the state):

MatrixClass getOperatorMatrix(size_t nrQubits, size_t qubit = 0, size_t controllingQubit1 = 0, size_t controllingQubit2 = 0) const override
{
    assert(qubit != controllingQubit1 && controllingQubit1 != controllingQubit2);

    const size_t nrBasisStates = 1ULL << nrQubits;
    MatrixClass extOperatorMat = MatrixClass::Zero(nrBasisStates, nrBasisStates);

    const size_t qubitBit = 1ULL << qubit;
    const size_t qubitBit2 = 1ULL << controllingQubit1;
    const size_t ctrlQubitBit = 1ULL << controllingQubit2;
    const size_t mask = qubitBit | qubitBit2 | ctrlQubitBit;

    // computing the tensor product between the gate matrix and identity operators for the other qubits

    for (size_t i = 0; i < nrBasisStates; ++i)
    {
        const size_t ind1 = i | mask;
        for (size_t j = 0; j < nrBasisStates; ++j)
            if (ind1 == (j | mask)) // the delta 'function'
                extOperatorMat(i, j) = BaseClass::operatorMat((i & ctrlQubitBit ? 4 : 0) | (i & qubitBit2 ? 2 : 0) | (i & qubitBit ? 1 : 0), (j & ctrlQubitBit ? 4 : 0) | (j & qubitBit2 ? 2 : 0) | (j & qubitBit ? 1 : 0));
    }

    return extOperatorMat;
}

MatrixClass getOperatorMatrix(size_t nrQubits, size_t qubit = 0, size_t controllingQubit1 = 0, size_t controllingQubit2 = 0) const override

{

assert(qubit != controllingQubit1 && controllingQubit1 != controllingQubit2);

const size_t nrBasisStates = 1ULL << nrQubits;

MatrixClass extOperatorMat = MatrixClass::Zero(nrBasisStates, nrBasisStates);

const size_t qubitBit = 1ULL << qubit;

const size_t qubitBit2 = 1ULL << controllingQubit1;

const size_t ctrlQubitBit = 1ULL << controllingQubit2;

const size_t mask = qubitBit | qubitBit2 | ctrlQubitBit;

// computing the tensor product between the gate matrix and identity operators for the other qubits

for (size_t i = 0; i < nrBasisStates; ++i)

{

const size_t ind1 = i | mask;

for (size_t j = 0; j < nrBasisStates; ++j)

if (ind1 == (j | mask)) // the delta 'function'

extOperatorMat(i, j) = BaseClass::operatorMat((i & ctrlQubitBit ? 4 : 0) | (i & qubitBit2 ? 2 : 0) | (i & qubitBit ? 1 : 0), (j & ctrlQubitBit ? 4 : 0) | (j & qubitBit2 ? 2 : 0) | (j & qubitBit ? 1 : 0));

}

return extOperatorMat;

}

For two and three qubits gates it’s more complex, but the idea is the same.
If you look into QubitRegister.h, QubitRegister::ApplyGate shows how easy it was using the matrix:

registerStorage = gate.getOperatorMatrix(NrQubits, qubit, controllingQubit1, controllingQubit2) * registerStorage;

1	registerStorage = gate.getOperatorMatrix(NrQubits, qubit, controllingQubit1, controllingQubit2) * registerStorage;

Yes, it was this easy. By the way, the RegisterStorage is quite easy. It could have had only the registerStorage member and NrQubits and the random stuff in there, but for speeding it up and for adding compute/uncompute ability, I added some more things in there.
Anyway, the optimized ApplyGate is more complex also due of my attempts to use open mp to speed up computations and also because it needs special handling for the different number of qubits gates. I even tried to convince the compiler to optimize better for diagonal and antidiagonal gates (qiskit aer for example goes much further, a diagonal gate might have 1 values, and it is foolish to multiply a lot of values with 1, isn’t it?). Here is how the implementation looks like for the one qubit gates (see QubitRegisterCalculator class):

static inline void ApplyOneQubitGate(const GateClass& gate, const VectorClass& registerStorage, VectorClass& resultsStorage, const MatrixClass& gateMatrix, const size_t qubitBit, const size_t NrBasisStates)
{
    const size_t notQubitBit = ~qubitBit;

    if (gate.isDiagonal())
    {
        for (size_t state = 0; state < NrBasisStates; ++state)
            resultsStorage(state) = state & qubitBit ? gateMatrix(1, 1) * registerStorage(state | qubitBit) : gateMatrix(0, 0) * registerStorage(state & notQubitBit);
    }
    else if (gate.isAntidiagonal())
    {
        for (size_t state = 0; state < NrBasisStates; ++state)
            resultsStorage(state) = state & qubitBit ? gateMatrix(1, 0) * registerStorage(state & notQubitBit) : gateMatrix(0, 1) * registerStorage(state | qubitBit);
    }
    else
    {
        for (size_t state = 0; state < NrBasisStates; ++state)
        {
            const size_t row = state & qubitBit ? 1 : 0;

            resultsStorage(state) = gateMatrix(row, 0) * registerStorage(state & notQubitBit) +
                gateMatrix(row, 1) * registerStorage(state | qubitBit);
        }
    }
}

static inline void ApplyOneQubitGate(const GateClass& gate, const VectorClass& registerStorage, VectorClass& resultsStorage, const MatrixClass& gateMatrix, const size_t qubitBit, const size_t NrBasisStates)

{

const size_t notQubitBit = ~qubitBit;

if (gate.isDiagonal())

{

for (size_t state = 0; state < NrBasisStates; ++state)

resultsStorage(state) = state & qubitBit ? gateMatrix(1, 1) * registerStorage(state | qubitBit) : gateMatrix(0, 0) * registerStorage(state & notQubitBit);

}

else if (gate.isAntidiagonal())

{

for (size_t state = 0; state < NrBasisStates; ++state)

resultsStorage(state) = state & qubitBit ? gateMatrix(1, 0) * registerStorage(state & notQubitBit) : gateMatrix(0, 1) * registerStorage(state | qubitBit);

}

else

{

for (size_t state = 0; state < NrBasisStates; ++state)

{

const size_t row = state & qubitBit ? 1 : 0;

resultsStorage(state) = gateMatrix(row, 0) * registerStorage(state & notQubitBit) +

gateMatrix(row, 1) * registerStorage(state | qubitBit);

}

This just takes advantage of the structure of that big matrix obtained by tensor product. How fast can it be? I tried it against qiskit aer with randomly generated circuits. On my computer, on windows, it beats qiskit aer (default setting, running on cpu, ‘statevector’, no config change) up to 18 qubits or so. For a small number of qubits it is much faster (as in several times faster). qiskit aer has optimisation not only with open mp, but also with AVX2 (Eigen has it too, but the operations done for applying a gate there are none)… and more, it has much better optimizations focused on particular kind of gates. There is more, it accumulates the gates without doing the computations and tries some optimizations when they need to be applied (before a measurement, for example), by joining and even eliminating some gates. I tried QCSim on up to 28 qubits. There QCSim can be even several times slower. But those are only some constants, one cannot improve more than O(N) in general (where N again is the number of basis states, 2^nrQubits or 1ULL << nrQubits as it’s computed in the project).

By the way, I had to use qiskit aer in a closed source project, that’s how I know, I had to look briefly even into its source code. Since I had some issues with compiling it on windows, I even did some fixes and became a contributor. Hopefully I’ll be of help there in the future, too.

Since I mentioned measurements, here is the implementation for the simpler one, a measurement of all qubits (for a single qubit or several ones cases I’ll let you check the code, it’s a little bit more complex):

size_t MeasureAll()
{
    const double prob = 1. - uniformZeroOne(rng); // this excludes 0 as probabiliy 
    double accum = 0;
    size_t state = NrBasisStates - 1;

    for (size_t i = 0; i < NrBasisStates; ++i)
    {
        accum += std::norm(registerStorage(i));
        if (prob <= accum)
        {
            state = i;
            break;
        }
    }

    setToBasisState(state); // collapse

    return state;
}

size_t MeasureAll()

{

const double prob = 1. - uniformZeroOne(rng); // this excludes 0 as probabiliy

double accum = 0;

size_t state = NrBasisStates - 1;

for (size_t i = 0; i < NrBasisStates; ++i)

{

accum += std::norm(registerStorage(i));

if (prob <= accum)

{

state = i;

break;

}

setToBasisState(state); // collapse

return state;

}

Conclusion

This merely scratched the surface. I presented briefly only the simulator parts of the project. There are plenty of algorithms implemented in there… luckily they can be investigated separately. For starting points, look into the cpp ‘Tests’ files, those show how the algorithms classes can be used, then look into the corresponding header files.

As usual, if you find issues in the project, please let me know.
Also it would be nice to give the project a star if you like it, since it might bring more exposure.

Later edit: Since I mentioned a couple of very good articles with tutorials, I added a repository for them on GitHub⁸

CodeProject

QCSim The QCSim repository on GitHub ↩
Quantum Computation and Quantum Information The Book by Michael Nielsen and Isaac Chuang ↩
Introduction to Coding Quantum Algorithms: A Tutorial Series Using Qiskit by Daniel Koch, Laura Wessing, Paul M. Alsing ↩
Fundamentals In Quantum Algorithms: A Tutorial Series Using Qiskit Continued by Daniel Koch, Saahil Patel, Laura Wessing, Paul M. Alsing ↩
Quantum Algorithm Implementations for Beginners by Abhijith J. at al ↩
Undergraduate computational physics projects on quantum computing by D. Candela ↩
Eigen The matrix library ↩
GitHub repository with the code for the tutorials articles mentioned above. ↩

Introduction

Algorithms

Theory

The code

Conclusion

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply

Introduction

Algorithms

Theory

The code

Conclusion

Share this:

Like this:

Related

Leave a Reply (for privacy concerns, please read the Privacy Policy Page accessible from the top menu)Cancel reply