量子連想メモリ@量子コンピュータ

量子コンピュータアルゴリズムである、http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.2291[quant-ph/9807053] Quantum Associative Memory (1998年)について。

「量子連想メモリ」という仰々しい名前がついていますが、{0,3,6,9,12,15}といったデータの中から、{6,7}を見つける問題です。この場合6がデータの中に入っているので6が見つかって欲しいです。

ソースコードは下の方にあります。QCLで実装しているので、量子コンピュータの基礎:振幅の初期化 - yukobaのブログを先にごらんください。アルゴリズムは、その「振幅の初期化 (1998年)」+「Groverのアルゴリズム (1996年)」です。

qcl -i quam.qcl を行い、testQuam() でテストコードを実行すると、

: SPECTRUM x: <0,1,2,3>
0.0026042 |0>, 0.0026042 |1>, 0.0026042 |2>, 0.0026042 |3>, 
0.0026042 |4>, 0.0026042 |5>, 0.9401 |6>, 0.023438 |7>, 
0.0026042 |8>, 0.0026042 |9>, 0.0026042 |10>, 0.0026042 |11>, 
0.0026042 |12>, 0.0026042 |13>, 0.0026042 |14>, 0.0026042 |15>
: found = 6 
[0/32] -1 |6>

のように結果が返ってきます。94%の確率で |6>、2.3%の確率で |7>、残りの確率でそれ以外が返ってきます。確率を100%に近づけるには、これを複数回反復すればOKです。

ソースコード

ファイル名は、quam.qcl

// Quantum Associative Memory (QuAM)

include "initAmpDist.qcl";
include "arcsin.qcl";

qufunct query(qureg x, quvoid f, int n) {
    int i;
    for i=0 to #x-1 {    // x -> NOT (x XOR n)
        if not bit(n,i) { Not(x[i]); }
    }
    CNot(f,x);            // flip f if x=1111..
    for i=0 to #x-1 {    // x <- NOT (x XOR n)
        if not bit(n,i) { !Not(x[i]); }
    }
}

operator flipOne(qureg q, quvoid f, int n) {
    query(q, f, n);
    CPhase(pi, f); // 符号反転 for f = 1
    !query(q, f, n);
}

operator diffuse(qureg q) {
    H(q);                        // Hadamard Transform
    Not(q);                        // Invert q
    CPhase(pi,q);                // Rotate if q=1111..
    !Not(q);                    // undo inversion
    !H(q);                        // undo Hadamard Transform
}

// Storing patterns in a quantum system
operator store(qureg x, int vector patterns) {
    int i;
    int vector values[#patterns];
    for i = 0 to #patterns - 1 {
        values[i] = 1;
    }
    initAmp(x, patterns, values);
}

// 論文の4.2
int loopCount(int xBitSize, int patternSize, real r1) {
    real N; real p; real r0;
    real a; real b; real k; real l;
    N = 2^xBitSize;
    p = patternSize;
    r0 = patternSize - r1;
    a = 2.0 * (p - 2.0 * r1) / N;
    b = 4.0 * (p + r0) / N;
    k = 4.0 * a - a * b + r1 / (r0 + r1);
    l = -a * b + 2.0 * a * (N + p - r0 - 2.0 * r1) / (N - r0 - r1) - (p - r1) / (N - r0 -r1);
    return floor((pi / 2.0 - arctan(k / l * sqrt((r0 + r1) / (N - r0 - r1)))) / 
        arccos(1.0 - 2.0 * (r0 + r1) / N));
}

// The main algorithm
operator modifiedGrover(int vector patterns, int vector searches, quvoid x) {
    int i; int j; int T;
    qureg f[1];
    
    // store
    store(x, patterns);

    // I_{\tau}
    for j = 0 to #searches - 1 {
        flipOne(x, f, searches[j]);
    }
    // G
    diffuse(x);
    // I_p
    for j = 0 to #patterns - 1 {
        flipOne(x, f, patterns[j]);
    }
    // G
    diffuse(x);

    T = loopCount(#x, #patterns, #searches);
    for i = 1 to T { 
        // I_{\tau}
        for j = 0 to #searches - 1 {
            flipOne(x, f, searches[j]);
        }
        // G
        diffuse(x);
    }
}

procedure quam(int patternBitSize, int vector patterns, int vector searches) {
    int found;
    qureg x[patternBitSize];

    reset;
    modifiedGrover(patterns, searches, x);
    dump x;
    
    measure x, found;
    print "found =", found;
}

procedure testQuam() {
    const patternBitSize = 4;
    int vector patterns = vector(0,3,6,9,12,15);
    int vector searches = vector(6,7);
    quam(patternBitSize, patterns, searches);
}

ファイル名は、arcsin.qcl

real arcsin(real x) {
    return Re(-I * log(I * x + sqrt(1 - x * x)));
}

real arccos(real x) {
    return Re(-I * log(x + I * sqrt(1 - x * x)));
}

real arctan(real x) {
    return Re(I / 2 * (log(1 - I * x) - log(1 + I * x)));
}

initAmpDist.qcl量子コンピュータの基礎:振幅の初期化 - yukobaのブログの通りです。