量子コンピュータの基礎：振幅の初期化

本記事は、量子コンピュータの初期状態の作り方についてです。1998年に発表された手法です。

お薦めの本

本題の前に、お薦めの本。量子コンピュータの書籍は色々出ていますが、下記の本がわかりやすかったです。原著、2007年。

マーミン 量子コンピュータ科学の基礎

ちなみに、素因数分解のアルゴリズムを開発し、量子コンピュータの火をつけた、ShorさんのMITの授業では、下記の本(原著、2000年)が使われています。

量子コンピュータと量子通信〈1〉量子力学とコンピュータ科学 (量子コンピュータと量子通信 1)

どの本を買っても、日本語で出ている本は、大旨、1996年までに発見されたアルゴリズムしか扱っていません。1997年以降のアルゴリズムは論文などを読む必要があります。本記事の内容(1998年)も同様です。

状態

2量子ビットの場合、0,1,2,3と4つの状態が作られます。を観測すると、振幅の2乗が出現確率なので、この状態を「観測」すると、50%の確率で|1>、50%の確率で|3>が得られます。

問題は、この状態をどうやって人為的に作り出すかというお話しです。このレベルですら、論文になるほどの内容…。論文は、http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.40.6693、[quant-ph/9807054] Initializing the Amplitude Distribution of a Quantum Stateです。アルゴリズムの説明を書こうかと思ったんですが、ややこしすぎるので省略。論文を読んでください。このアルゴリズムは、希望した状態のところに等確率を割り振るだけです。2:1:1とかはできません。ただし、確率を保持したまま、位相をずらすことはできます。2:1:1の作り方は、量子コンピュータで自由に初期状態を作る方法 - yukobaのブログをご覧ください。

その代わりにソースコードを載せます。プログラミング言語はQCLです。QCL - A Programming Language for Quantum Computers。これ、面白いことに、古典コードと量子コンピュータのコードを混在することができます。将来、両タイプのチップが載った時に、自動的に処理を割り振るように実装されることを想定して作っています。ちなみに、QCLはUbuntu 8.04では動いたんですが、cygwin での動かし方がわかりませんでした。Ubuntu 8.04 の場合、少なくとも、libplot-devパッケージを事前にインストールする必要があります。

論文には、位相のずらし方は、N段階でずらす方法も書いてあるのですが、最初の方に書いてある、1 または -1 の2段階のずらし方のみの実装です。qcl -i initAmpDist.qcl で起動し、testInitAmp() でテストコードを実行できます。

: SPECTRUM x: <0,1,2>
0.33333 |0>, 0.33333 |2>, 0.33333 |5>
[0/32] 0.57735 |0> - 0.57735 |2> - 0.57735 |5>

という結果が返ってきます。|0>, |2>, |5> にそれぞれ、1/3の確率をふり、|2>と|5>の位相は-1倍しています。

量子コンピュータで自由に初期状態を作る方法 - yukobaのブログへと続く。

ソースコード

では、ソースコード。ファイル名は、initAmpDist.qcl。

// Initializing amplitude distribution (振幅分布の初期化)
// 論文の上位ビットと下位ビットの順番が一般の順序と逆であることに注意！
// 特にcのビット順が問題。

// 論文の3〜6行目
qufunct initAmpFlip(qureg x, qureg c, int pattern, int prevPattern) {
  int j;
  // 論文の5,6行目はF^0であり、CNotはF^1なので、ビット反転。
  Not(c[0]);
  // 論文の3〜5行目
  for j = 0 to #x - 1 {
    if bit(pattern, j) xor bit(prevPattern, j) { CNot(x[j], c[0]); }
  }
  // 論文の6行目
  CNot(c[1], c[0]);
  // 反転したc[0]を元に戻す
  Not(c[0]);
}

// 論文の7行目。s = {1, -1}。
operator initAmpStateGen(qureg c, int s, real p) {
  real p1; real p2;
  // p1^2 + p2^2 = 1 より行列はユニタリ
  p1 = sqrt((p - 1) / p);
  p2 = s / sqrt(p);

  Matrix4x4(
    1, 0, 0, 0,
    0, 1, 0, 0,
    0, 0, p1, -p2,
    0, 0, p2, p1,
    c);
}

// c1 == i1 かつ c2 == i2 の時に、a を反転する。i1, i2 = {0, 1}。#a, #c1, #c2 = 1。
qufunct fredkin(qureg a, qureg c1, boolean i1, qureg c2, boolean i2) {
  if not i1 { Not(c1); }
  if not i2 { Not(c2); }
  CNot(a, c1 & c2);
  if not i2 { Not(c2); }
  if not i1 { Not(c1); }
}

// 論文の8〜10行目
// x が z と等しければ、g[#x - 2] のフラグを立てる
qufunct initAmpAnd(qureg x, quvoid g, int z) {
  int k;
  fredkin(g[0], x[0], bit(z, 0), x[1], bit(z, 1));
  for k = 2 to #x - 1 {
    fredkin(g[k - 1], x[k], bit(z, k), g[k - 2], true);
  }
}

// 振幅分布の初期化。
operator initAmp(quvoid x, int vector patterns, int vector values) {
  int i; int prevPattern;
  qureg g[#x - 1];
  qureg c[2];
  
  prevPattern = 0;
  for i = 0 to #patterns - 1 {
    // 論文の3〜6行目
    initAmpFlip(x, c, patterns[i], prevPattern);
    prevPattern = patterns[i];
    
    // S^p。論文の7行目。
    // initAmpStateGen(c, 1, #patterns - i);
    initAmpStateGen(c, values[i], #patterns - i);
    
    // save。パターンに等しいやつのc[1]を0に戻す。
    // 論文の8〜10行目(AND)
    initAmpAnd(x, g, patterns[i]);
    // 論文の11行目
    CNot(c[1], g[#x - 2]);
    // 論文の12〜14行目(AND逆戻し）
    !initAmpAnd(x, g, patterns[i]);
  }
  // 論文の15行目
  Not(c[0]);
}

procedure testInitAmp() {
  const patternBitSize = 3;
  int vector patterns = vector(0, 2, 5);
  int vector values = vector(1, -1, -1);
  qureg x[patternBitSize];
  initAmp(x, patterns, values);
  dump x;
}

http://d.hatena.ne.jp/yukoba/20090929/p1へと続く。