学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月

11-NA14量子 i.i.d. 状態のシミュレーション
坂下達哉（電気通信大学）

概要
量子情報理論において，量子 i.i.d. 状態（n 次のテンソル積で表現される量子状態）に関する極限式は重要である．本研究では， 二つの 2 × 2 の密度行列から成る量子 i.i.d. 状態の仮説検定問題における誤り確率を扱う．これを直接的な方法で計算しようとすると，2n × 2n という非常に大きなサイズの行列の固有値問題を解く必要が生じる．しかし，テンソル積の既約分解という数学的概念にもとづいた計算手法を用いると，最大サイズ (n + 1) × (n + 1) までのいくつかの行列の固有値問題を解くことによって同じ計算が可能になる．本稿では，この計算方法に関して本年度取り組んできた並列実装上の工夫について紹介する．
した [10, 15]．本研究の目的は，既知の問題を解析し，未解決の問題の極限式予測に活用でき

1 研究の目的と意義るようにすることである．例えば，文献 [10] では，中心極限定理に関連した新しい予想を数値量子情報理論を理解する上で，量子 i.i.d. 状的に検証している．

態に関する極限式は重要である．これまで，古典情報理論において漸近性を扱うには数値的手

本年度は，このような新しい予想への適用に法が活用されていたが，量子情報理論において備えて，既知の仮説検定問題を題材として並列は実現可能とも有効であるとも考えられていな実装上の工夫を行うことに専念した．以下の副かった．

節では，その工夫について述べるための準備として，研究の概要を紹介し，これまで得られたテンソル積の既約分解という手法は林 [3] に成果について述べる．

より初めて量子 i.i.d. 状態に適用され，長岡 [8]により数値計算アルゴリズムが提案された．これにより，数値計算が現実的な時間で行えるよ

1.1 量子仮説検定うになった．このアルゴリズムは，長岡の指導のもとで，柿崎 [12]，堂嶋 [16] により研究さ

本節では量子仮説検定について必要最小限のれた．

事柄を述べる．

これらの研究を受けて，我々は，既知の量子ま ず ，い く つ か の 物 理 的 な 概 念 に 数 学 的

仮説検定問題を題材として，大規模で高精度のな 定 義 を 与 え る ．量 子 状 態 は ，ト レ ー ス が数値シミュレーションが可能であることを実証

1 の半正定値行列（密度行列と呼ばれる）ρ
1

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月

で 表 さ れ る ．測 定 は ，半 正 定 値 行 列 の 集 合と，A のスペクトル分解∑

M = {M (x)}x∈X で，x∈X M (x) = I（単位A=行列）を満たすもので表される．これを POVM

∑λi xi x∗i(4)i(Positive Operator-Valued Measure) と呼ぶ．

が得られる．ここで，xi は縦ベクトルで x∗i はここで，添字 x は測定値を表す．量子状態 ρ の横ベクトルであるから，積 xi x∗i は A と同じサもとで測定 M を行ったとき，測定値 x が得らイズの正方行列である．これに対し，れる確率は{A > 0} :=PρM (x) := Tr[ρM (x)]∑xi x∗i(5)λi >0(1)と定義する．{A > 0} は A と同じサイズのエルミート行列であり，A の正固有値に対応するで表される．

固有空間の直和への射影を表す．

量子 i.i.d. (independent and identically dis-式 (3) は量子 Neyman-Pearson 検定 {ρ⊗n −tributed) 状態は，n 次のテンソル積

ena σ ⊗n > 0} の二つの誤り確率のうちの一つ

ρ⊗n = ρ ⊗ · · · ⊗ ρ| {z }であり，σ ⊗n が真のときに ρ⊗n を誤って真とn判断する確率を意味する [8]．また，で表される状態のことである．ここで，行列′1log βn (a),nr(a) := lim rn (a)rn (a) := −′A : d × d サイズと B : d × d サイズのテンソ
ル積を dd′ × dd′ サイズ の行列a11 B ..
A ⊗ B :=  .
ad1 B···
..
.
···

n→∞


a1d B.. . (6)(7)とすると，量子 Hoeﬀding の定理 [4, 6] より，任意の a > −D(σ∥ρ) に対して，(2)r(a) = max{θa − ψ(θ)}

add Bθ∈R(8)と定義する．ただし，n はテンソル次数とよばここで，ψ(θ) = log Tr[ρθ σ 1−θ ]れる．量子 i.i.d. 状態は，物理的には独立に用が成り立つ．式 (8) の右辺は Newton 法などを意された n 個の状態を意味しており，古典確率用いることで求められ，計算する行列は 2 × 2論における独立同一分布に対応する．

サイズなので計算時間はかからない．

次に，本研究で扱う量子仮説検定問題につい量子情報理論では，式 (7) のような量子 i.i.d.て述べる．いま，2 × 2 サイズの密度行列 ρ, σ状態に関する n → ∞ における極限式が至るで表される量子状態が与えられているとし，実ところに現れる．冒頭でも述べたが，本研究の数 a に対して，次の量目的は，このような極限式の予想を行うために

βn (a) := Tr[σ ⊗n {ρ⊗n − ena σ ⊗n > 0}] (3)数値シミュレーションを活用することである．

を考える．ここで，任意のエルミート行列 A に限式 (8) が得られている量子 Hoeﬀding の定理

対し，行列 {A > 0} を以下のように定める．A

は，プログラムの試作に当たって格好の題材での固有値と固有ベクトルから成る正規直交基底ある．そこで，以降では，量子 Hoeﬀding の定をそれぞれ {λi }，{xi }（縦ベクトル）とおく

理に的を絞って解説を行う．

そのための並列化実装を工夫する上で，既に極2

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月

αk,ij (i, j ∈ {0, 1, . . . , r}) は，まず1.2 既約成分の計算アルゴリズム
r−j⊗n

量子 i.i.d. 状態の密度行列 ρ(a11 y+a21 z)のサイズは′αk,ijy r−i z i(12)n′という多項式展開から係数 αk,ijを求め，これは大変である．そこで，本節ではテンソル積のを用いて既約分解という方法を導入する．

√任意の 2 × 2 行列

[A=a11a21a12a22αk,ij = (det A)k](9)OA1(13)である．よって，テンソル積の加算・減算は既Ak(10)約成分ごとに行える．また，後でみるようにトレースをとるので，この基底変換は存在するこ



=



′
αk,ij

いる基底変換が A の行列成分によらないこと

k=0

A0
r Ci

ここで，肝要なことは式 (11) の既約分解に用⌊n/2⌋ mk⊕ ⊕r Cjと計算できる．

に対し適当な基底変換を施すことで，A⊗n はO

このように，テンソル積の既約分解を用いる


A1..とは重要であるが具体的に求める必要はない．

.

ことで，指数サイズであるテンソル積の行列成分の計算が多項式オーダーで行える．

1.3 既約分解を用いた誤り確率の計算

A⌊ n ⌋2

本節では前節で述べたテンソル積の既約分解(11)を誤り確率の計算に適用する．
式 (11) の既約分解を用いることによって，という形のブロック対角行列で表される [2]．

誤り確率 βn (a) は次のように表される:

各 Ak はテンソル積 A⊗n の既約成分とよばれ，そのサイズは dim Ak = n + 1 − 2k で，βn (a) = Tr[σ ⊗n {ρ⊗n − ena σ ⊗n > 0}] (14)

重複度 mk := n Ck − n Ck−1 回だけ現れる．

⌊n/2⌋=（ここで，便宜的に n C−1 := 0 とする．）つ∑mk βn,k (a)(15)βn,k (a) := Tr[σk {ρk − ena σk > 0}]

(16)

k=0

ま り ，式 (11) は ま ず A0 を m0 = 1 個 だ
ここで，け 対 角 に 並 べ ，次 に A1 を m1 個 だ け 対 角に並べるという操作を繰り返して得られる∑(a12 y+a22 z) =r∑
i=02 × 2 になるため，この行列についての計算
nj

nk (n Ck − n Ck−1 )(n + 1 − 2k) = 2 サイズの正方行列である．ここで，便宜的に n C−1 := 0であり，式 (16) の {ρk } と {σk } はそれぞれとしている．

ρ⊗n と σ ⊗n の既約成分である．

r := n − 2k とおくとき，Ak の (i, j) 成分

図 1 で式 (16),(15) の計算手順を与える．

3

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月（ステップ）1. 各 k について，以下を行う．

(a) 式 (12),(13) を用いて {ρk } と {σk } を計算する．
(b) 行列 ρk − ena σk の固有値 {λn,k,i }i と正規直交化された固有ベクトル {vn,k,i }iを求める．

(c) βn,k (a) =∑∗vn,k,iσk vn,k,i を求める．

λn,k,i >0⌊n/2⌋∑（ステップ）2. βn (a) =mk βn,k (a) を求める．

k=0

図1βn (a) の計算手順

ここで，図 1 の手順のステップ 1 において，各 k に対する処理は独立なので並列化を施し1.4 計算機環境ている．また，複数の a について βn (a) を求める場合，図 1 のステップ 1.(a) において既約成使用した計算機環境は以下の通りである．

分 {ρk }, {σk } は全ての a について共通なので，• 計算機 : スーパーコンピュータ 東京大

一度だけ計算すればよいことに留意する．以上の二点についての詳細は 1.5 節で述べる．

学 HA8000 (AMD Opteron processor

2 × 2 サイズの２つの密度行列については，8356)

実対称行列に限って良いことが知られているの• コンパイラ : Intel C++ compiler

で，実際の計算では ρ, σ は実対称行列のみを• C++ テンプレート線形代数ライブラリ用いる．固有値分解にも実対称行列向けのルー: Eigen3[1]チンを用いる．固有値ルーチンを選定する際に• 多倍長演算ライブラリ : Exflib[11]は，以下に事項に注意を要する．

• 固有値分解に用いるライブラリ : LAPACK dsyev 関数 (倍精度・QR 法の場

1. 全固有値，全固有ベクトルが必要である．

合) ，Eigen3 の固有値分解関数を改変

2. 行列 ρ⊗n , σ ⊗n は指数的にばらついた固(多倍長・QR 法の場合)有値を持つ．したがって，実対称行列

ρk − ena σk には絶対値の小さい正負の工夫した点としては，Eigen3 に多倍長型を固有値があるが，固有値の正負が数値誤追加し，多倍長型を成分に持つ行列を扱える差によって変わることがないように正確ようにしたことがあげられる．これにより，倍に求める必要がある．

精度に限らず多倍長も MATLAB などと同様に行列演算を記述できる．また，倍精度・多倍

このような場合，固有値解法には Householder

長の切り替えも，行列クラス名の置換のみで済QR 法が適していることがわかっているためむ．なお，現在，上述の Eigen3 の固有値分解[14]，これを用いる．

関数の性能は不十分なため，自前で改良し並列化を施すことも本研究の課題である．

4

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月めに，MPI_Gather 関数を用いることにする．
ここで，計算結果の集約と総和の両方を行える

1.5 これまでの並列化方法

MPI_Reduce 関数を用いないのは，mk βn,k (a)本節では，昨年度までの並列化方法についても考察の対象にするので，これもデータとして記述する．次節で，この実装の問題点と今回新出力したいからである．

しく考案した並列化法を述べる．

最 後 に ，ル ー ト プ ロ セ ス に お い て ，誤

図 1 の計算手順において，ステップ 1 の異なり 確 率 βn (a) を 用 い て 誤 り 指 数 rn (a) =− log βn (a)/n を計算する．

る k についての処理は互いに独立であるため，自明な並列化が施せることを述べた．

例として，テンソル次数 n = 100 の場合を並列化方法としては，各プロセス（スレッド）とりあげる．このとき，既約成分の番号 k はに複数の k についての処理を担当させる必要

0, 1, · · · , ⌊n/2⌋ = 50 で，対応する既約成分のがある．その際，各既約成分の格納には膨大なサイズ dim Ak は 101, 99, . . . , 3, 1 となる．プ

メモリが必要で，データの送受信は各プロセスロセス数が 4 の場合について，ρ⊗n と σ ⊗n の（スレッド）で計算した mk βn,k (a) と固有値分既約成分の格納方式と計算結果の集約方法を図解の計算時間等の付加的な情報だけである．こ2，図 3 に示す．ここで，ルートプロセスはプのような処理には OpenMP よりも MPI が向ロセス 0 としている．

いているため，本研究では MPI を用いることにする．

2 当拠点公募型共同研究として各プロセスのメモリ使用量・計算量をできる実施した意義

だけ均一にするため，担当する既約成分は巡回的に割り当て，この割り当てはすべての a に対して固定しておく．また，既約成分 {ρk }, {σk }(1) 共同研究を実施した大学名：東京大学

は一度だけ計算してメモリに格納しておき使い(2) 共同研究分野：超大規模数値計算系応用回すことにする．以下で，一つの a に対して誤分野

り確率 βn (a) の具体的な計算手順を述べる．

(3) 当公募型共同研究ならではという事項な

まず，各プロセスで既約成分を計算し，多倍ど：長型（または倍精度型）の行列に収める．その数値計算・並列化の専門家との協業が行え際，既約成分のサイズは各々異なるので，行列たことは有意義であった．特に，固有値解クラスからなる配列を用いて管理することはで法における三重対角化の並列化方法についきない．そこで，「行列クラスへのポインタ」かてアドバイスをいただくことができた．まらなる配列を用いて既約成分を管理する．

た，スーパーコンピュータ HA8000 は In-

次に，各プロセスで mk βn,k (a) を計算し，tel 互換の CPU を採用しているので，手元あらかじめ決めておいた一つのプロセスのパソコンで試したプログラムをそのまま（ルートプロセスと呼ぶ）に集約した後，ルー使用できるという利点がある．特に，C++ト プ ロ セ ス で ス テ ッ プ 2 の 総 和 βn (a) =のテンプレートライブラリ Eigen3 は対応

k mk βn,k (a) を計算する．この集約を行うた

していない CPU もあるが，HA8000 では∑

5

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月
䠌䠓䠍䠌䠍䠌䠍ḟඖ䠔䠓ḟඖ䠔䠍ḟඖ䠍䠐䠍䠍䠕䠕ḟඖ䠕䠏ḟඖ䠓䠕
ḟඖ䠎䠔

͙䠑䠌

⾜ิ䜽䝷䝇䜈䛾䝫䜲䞁䝍䛾㓄ิ᪤⣙ᡂศ䛾ᐤ୚䛸䛧䛶ฟຊ᪤⣙ᡂศ䛾␒ྕ᪤⣙ᡂศ䛤䛸䛻ㄗ䜚☜⋡䛾ィ⟬ DW/ͺ'ĂƚŚĞƌ⾜ิ䜽䝷䝇䜈䛾䝫䜲䞁䝍䛾㓄ิ䠍ḟඖ͙䠋䠌䠓䠍䠌͙䠐䠕䠍䠐䠍䠍͙䠑䠌䠎䠑䠔͙䠋䠏䠒䠕͙䠐䠔䠌⾜ิ䜽䝷䝇䜈䛾䝫䜲䞁䝍䛾㓄ิ

m= 0 1 2䠕䠍ḟඖ䠔䠑ḟඖ䠏䠒䠕͙䠔䠏ḟඖ䠔䠕ḟඖ䠔䠕ḟඖ䠐䠕䠏ḟඖ䠕䠓ḟඖ䠕䠑
ḟඖ

図2䠑͙

⧞䜚㏉䛧ᩘ
䠐䠔

⾜ิ䜽䝷䝇䜈䛾䝫䜲䞁䝍䛾㓄ิn/䝥䝻䝉䝇ᩘ䠑ḟඖ図3ρ⊗n と σ ⊗n の既約成分の格納方式

䠍͙䠑䠌

⥲࿴䜢䛸䜛
1− logn➨஧✀ㄗ䜚䝺䞊䝖

今までの計算結果の集約方法• 支配的な既約成分を精度よく求める．

問題なく使用できた．

• その他の既約成分で大きな数値誤差が起こらないようにする．

3 研究結果の詳細と当初計画の以上を実現するために，次のように実装およ達成状況

び MPI による並列化を変更した．

1. 支配的な既約成分だけを重点的に計算し，不3.1 研究結果の詳細について要な既約成分の計算を省略する．そのため，小さな n についてまず計算することで，支

本節では，平成 23 年度のプロジェクトにお配的な既約成分を予想する．

いて取り組んだ以下の課題について記述する．

2. これまでは，全ての既約成分についての計(課題 1) 並列化の実装方法についての工夫（中算が終わってから集団通信を用いて一箇所間報告までに実施した課題）に集める方式をとっていたため，プロセス(課題 2) 固有値計算を行う際の三重対角化の間の待ち時間が生じていた．そこで，非同OpenMP によるスレッド並列化 （最期通信を用いることにより，計算と通信の終報告までに実施した課題）オーバーラップを図る．

課題 1 は以下に述べる「支配的な既約成分」

3. 大きな既約成分ほど数値誤差が大きいことという知見に基づく．式 (15) により誤り確がわかっている．そのため，小さな既約成率 βn (a) は各既約成分 k についての計算値分から計算し大きな既約成分で数値誤差をmk βn,k (a) の総和として求められる．我々は，検出したら計算を打ち切る（図 5 参照）．

その総和において実質上，ごく一部の項のみが以上について実装は出来上がっているが，現寄与することを見出した [14]．これらの項に関段階において上の実装を施す意義は必ずしも明する既約成分を支配的な既約成分と呼ぶ（図 4

確ではない．そのため，この実装が生かせる問参照）．

題設定を提案して，その問題に適用した結果と支配的な既約成分を考慮すると，以下の点にともに対外発表に望みたい．

留意すれば効率の良い計算が行える．

6

学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月
0.0008n = 100
n = 110n = 120n = 130
n = 1400.00070.0006䠌

0.0005᪤⣙ᡂศ䛾␒ྕ

0.0004
䠍

䠐

䠍䠍

䠓

䠍䠌

͙

䠑䠌

䠌

䠍

䠐䠕

͙

ィ⟬䛾೵Ṇ䞉
⥅⥆䛾ᣦ♧ィ⟬⤖ᯝ䛾䝁䝢䞊䠎䠏͙

0.0003䠐䠓͙

0.0002

䠎䠑䠔

͙䠋

䠏

䠒

䠕͙䠐䠔

䠐䠔

䠐䠕䠑䠌⥲࿴䜢䛸䜛ㄗ䜚☜⋡

0.0001000.1

0.2

0.30.4

支配的な既約成分
既約成分
大のサイズ

− n1 logィ⟬⤖ᯝ
䛾㏦ಙ

0.5

ㄗ䜚ᣦᩘ

小

図 5 新しい並列化方法

図4

n を変えたときの既約成分ごとの計算値(a = −1.5)

課題 2 を行う動機について述べる．現時点
で は n が 小 さ い と き ，HA8000 で 使 用 可 能

3.2 当初計画の達成状況についてな 64 × 16 = 1024 コアが生かしきれていな
い．そこで，多倍長の Householder QR 法をライブラリの解読と書き換えに多大な時間をOpenMP により並列化を行う．
Householder QR 法は，Householder 法によ要し，課題 2 の実装の完成に至らなかった．課

る三重対角化と三重対角行列の固有値・固有ベ

題 2 の実装で起こっている問題を解決するためクトルを求める QR 法からなる．本研究では，に，テンプレートライブラリの内部構造についHouseholder 三重対角化の OpenMP による並てさらに理解を深め，ライブラリの作成者とも列化に取り組んだ．

連絡をとり解決策を考えたい．

まず，三重対角化のうち，rank2update の for
文を OpenMP により並列化することから始め4 今後の展望た．テンプレートライブラリ Eigen3 のソースコードを解読し三重対角化に必要な部分のみを取り出した．Eigen3 はどのような備え付けの5 研究成果リスト型（float 型，倍精度等）でも受け付けるようになっているが，簡単のためこれを倍精度に書き研究成果のリストは以下の通りである．

直した．倍精度の場合には，OpenMP による並列化は成功した．だが，多倍長の場合につい(1) 学術論文 なして OpenMP による並列化を行うと，２スレッ(2) 国際会議プロシーディングス なし（発表ド以上で明らかに誤った計算結果となる．この予定）

原因について，現在調査を行っている．

(3) 国際会議発表 なし(4) 国内会議発表 なし(5) その他（特許，プレス発表，著書等） なし
7学際大規模情報基盤共同利用・共同研究拠点 平成 23 年度共同研究最終報告書 2012 年 5 月

sical and quantum hypothesis testing forsimple hypotheses.

参考文献

IEEE Trans. Inf.

Theory, 53(2):534–549, 2007.
[1] 線 形 代 数 計 算 用テ ン

[9] M. Nussbaum and A. Szkola. The Cher-

Eigen2.

noﬀ lower bound for symmetric quan-

http://eigen.tuxfamily.org/index.

tum hypothesis testing. Ann. Statist.,

php?title=Main_Page.

37(2), 2009.

C++

プ レ ー ト ラ イ ブ ラ リ[2] R. Goodman and N. R. Wallach. Repre-[10] T. Sakashita and H. Nagaoka. A nu-sentations and Invariants of the Classi-merical study of hypothesis testing forcal Groups. Cambridge University Press,

quantum i.i.d. states. Asian conference

1998.

on quantum information science(AQIS),[3] M. Hayashi. Asymptotics of quantum

pp. 213–214, Aug. 2010.

relative entropy from a representation[11] 藤 原.

theoretical viewpoint. J. Phys. A: Math.

並列化対応多倍長ライブ

ラ リ Exflib.

Gen., 34(16):3413–3419, 2001.

http://www-an.acs.i.

Error exponent inkyoto-u.ac.jp/~fujiwara/exflib/.
[12] 柿崎. 量子的 i.i.d. 情報源に関する数値計

asymmetric quantum hypothesis test-算によるアプローチ. Master’s thesis, 電ing and its application to classical-

気通信大学情報システム学研究科, 2006.

quantum channel coding. Phys. Rev. A,[13] 坂下. 量子 i.i.d. 状態のシミュレーション[4] M. Hayashi.

76(6):062301–1,062301–4, 2007.
[5] K. M. R. Audenaert et al.とその理論的考察. スーパーコンピュー

Discrimi-ティングニュース, Vol.13 特集号 1, Mar.

nating states: The quantum Chernoﬀ

2011.

bound. Phys. Rev. Lett., 98(16):160501–[14] 坂下, 片桐, 長岡. 量子 i.i.d. 状態の仮説検

1,160501–4, 2007.

定に関する数値的手法とその誤差分析. 情[6] H. Nagaoka. The converse part of the報処理学会論文誌コンピューティングシ

theorem for quantum Hoeﬀding bound.

ステム (ACS), (35), Oct. 2011.

arXiv:quant-ph/0611289v1, 2006.

[15] 坂下, 長岡. 量子 i.i.d. 状態のシミュレー[7] H. Nagaoka. Some mathematical prob-ションとその理論的考察. 日本応用数理

lems related to quantum hypothesis test-

学会 2010 年度年会 講演予稿集, pp. 245–

ing. In General Theory of Information

246, 2010.

Transfer and Combinatorics, Vol. 4123

[16] 堂嶋, 片桐, 長岡. 量子 i.i.d. 状態におけof Lecture Notes in Computer Science,

る仮説検定の漸近特性に関する数値的ア

pp. 1100–1103. Springer Berlin / Heidel-プローチ. 第 29 回情報理論とその応用シ

berg, 2006.

ンポジウム予稿集 (SITA), pp. 767–770,[8] H. Nagaoka and M. Hayashi.

An2006.

information-spectrum approach to clas8