2022 年度科学研究費シンポジウム

大規模複雑データの理論と方法論〜新たな発展と関連分野への応用〜
科学研究費補助金 基盤研究 (A) 20H00576「大規模複雑データの理論と方法論の革新的展開」（研究代表者:青嶋誠），学術研究助成基金助成金 挑戦的研究 (萌芽) 22K19769「テンソル構造をもつ巨大データの統計的圧縮技術の開発」（研究代表者: 青嶋誠）によるシンポジウムを下記のように催しますので，ご案内申し上げます．

誠(筑波大学)

矢田和善(筑波大学)

中山優吾(京都大学)

青嶋記
日時：2022 年 11 月 4 日 (金)〜5 日 (土)場所：つくば国際会議場 中会議室 202 ( https://www.epochal.or.jp/ja/ )〒305-0032 茨城県つくば市竹園 2-20-3

プログラム
11 月 4 日 (金)13:00〜13:10

開会

13:10〜13:50

藤森洸

（信州大・経法学部）

後藤佑一（九州大・大学院数理学研究院）

劉 言
谷口（早稲田大・理工学術院）
正信（早稲田大・理工学術院）

高次元・定常時系列に対するスパース主成分分析

13:55〜14:35

赤間陽二（東北大・理学研究科）

Unbounded largest eigenvalue of large sample correlation matrices:
Asymptotics and applications
14:40〜15:20

永井勇（中京大・教養教育研究院）

GMANOVA モデルにおける新たな推定方法と解釈
15:35〜16:15

中川智之 （東京理科大・理工学部）

ロバストダイバージェンスを用いたベイズ推論について

16:20〜17:00

寺田

吉壱（大阪大・大学院基礎工学研究科）

山本

倫生（大阪大・大学院人間科学研究科）

代表点を用いた大規模クラスタリングの近似法とその性質

17:05〜17:45

竹内努（名古屋大・理学研究科）

Recent progress in the application of high-dimensional statistics toastrophysics and cosmology

11 月 5 日 (土)9:00〜9:30

岡崎彰良（電気通信大・大学院情報理工学研究科）

川野秀一（九州大・大学院数理学研究院）

凸クラスタリングによるマルチタスク学習

9:35〜10:05 書川 侑子 （電気通信大・大学院情報理工学研究科）
川野秀一（九州大・大学院数理学研究院）

馬蹄事前分布に基づく連結 lasso ロジスティック回帰モデリング

10:10〜10:50

仲北祥悟（東京大・大学院総合文化研究科）

凸最適化によるエルゴード的拡散過程のオンライン推定

10:55〜11:35

今泉允聡（東京大・大学院総合文化研究科）

非スパースな高次元漸近論の理論と応用

11:35〜13:00

13:00〜13:40

昼食

橋本真太郎（広島大・大学院先進理工系科学研究科）

羽村 靖之（京都大・大学院経済学研究科）

鬼塚 貴広（広島大・大学院先進理工系科学研究科）

菅澤 翔之助（東京大・空間情報科学研究センター）

Sparse Bayesian inference on gamma-distributed observation
13:45〜14:25

矢野
清

恵佑
智也

（統計数理研究所）（東京大・大学院情報理工学系研究科）

Minimum information dependence modeling for mixed domain data
14:30〜15:10

鈴木大慈（東京大・大学院情報理工学系研究科）

高次元データ学習における特徴学習の優位性

15:25〜16:05 鹿野 豊

（群馬大・大学院理工学府, Institute for Quantum Studies, Chapman University）

量子乱数プロトコルを用いた量子計算機の安定性評価

16:10〜16:50

福地純一郎（学習院大・経済学部）

Bootstrap for selecting a subset which contains all populations better than a standard
16:50〜17:00

閉会高次元・定常時系列に対するスパース主成分分析
藤森洸 1

後藤佑一 2

劉言 3, 4

谷口正信 3

信州大学 経法学部
九州大学 数理学研究院
3
早稲田大学 基幹理工学部
理工学術院総合研究所 重点研究領域「数理科学研究所」
1

2

4

1

早稲田大学

概要

主成分分析は多変量データの共分散構造の解析に際し，中心的な役割を果たす重要な手法である．近年は，データの次元 p がサンプル数 n と比較して大きい場合，すなわち高次元の設定における主成分分析の研究が盛んである．高次元の設定において，従来の主成分分析は精度が悪いことが知られている．特に，スパースな主成分ベクトルの推定には，より良い性質を持った推定量が，主に独立標本の設定において提案されてきた．本研究は，高次元の設定におけるスパースな第一主成分の推定問題を，定常時系列データに対して適用し，得られる推定量の誤差評価を行うことを目的としている．具体的には，Gaussian よりも裾の重いケースを含む高次元定常時系列に対して Lasso 型のスパース主成分分析手法を適用し，得られる推定量の漸近的な挙動について議論する．さらに，従来の主成分分析手法や独立標本に対する Lasso 型主成分分析との比較を行っていく．

2

主結果{Xt }t∈Z を確率空間 (Ω, F, P ) 上の Rp -値，平均 0 の定常過程とする．観測系列 X1 , . . . , Xn , n ∈N に対して，次の p × p 標本共分散行列及び，共分散行列を考える．
1∑
Xt Xt⊤ ,
n
n

Σ̂n =

Σ0 = E[Xt Xt⊤ ].

t=1

以下，Σ0 の最大固有値 ϕ2max に対応する，規格化された第一主成分ベクトル q 0 の推定問題
を考える．推定の対象となるパラメータは，β 0 = ϕmax q 0 として特徴づけることができる
が，これは，次の最適化問題の解として与えられる：1
β 0 = arg min ∥Σ0 − ββ ⊤ ∥2F ⇐⇒ Σ0 β 0 = ∥β 0 ∥22 β 0 ,β 4ただし，∥·∥F は行列のフロベニウスノルムである．ここで，添え字集合 S を S = {j : βj0 ̸= 0}で定義する．特に，β 0 は s0 -sparse, つまり，|S| = s0 であると仮定する．さて，β 0 の推定量として，次を定義する (van de Geer (2016))：Deﬁnition 2.1.

{

β̂n := arg min
β∈B

}
1
⊤ 2
∥Σ̂n − ββ ∥F + λ pen(β) ,4

B := {β : ∥β − β 0 ∥2 ≤ η},

ただし，λ ≥ 0 は tuning parameter, pen(·) は罰則関数，η > 0 は適当な定数である．特に，本研究では次の Lasso 型推定量に焦点を当てる．
{
}
1
1
⊤ 2
β̂n := arg min∥Σ̂n − ββ ∥F + λ1 ∥β∥1 ,
(2.1)
β∈B
4
{
}
1
β̂n0 := arg min∥Σ̂n − ββ ⊤ ∥2F + λ0 ∥β∥0 .
(2.2)
β∈B
4

Lasso 型推定量を用いて，最大固有値 ϕ2max 及び対応する規格化された固有ベクトル q 0の推定量を次のように構成することができる：
ϕ̂2max := ∥β̂n1 ∥22 ,q̂n :=

β̂n1∥β̂n1 ∥2.

これらの推定量について，n → ∞ のとき，mixing 条件や，スパイク性に関する条件をはじめとする適当な正則条件の下で，α-mixing Gaussian process に対する推定量について次が成立することが示される：
(
)
√
)
(
log
p
log p
3/2
1
0
1
0 2
3 2
∥β̂n − β ∥1 = Op s0 ϕmax,
, ∥β̂n − β ∥2 = Op s0 ϕmaxnn(|ϕ̂max − ϕmax | = Op√
3/2
s0 ϕmax

log p
n

)

(
,

∥q̂n − q ∥2 = Op
0

√
3/2
s0

log p
n

)
.

講演内では，Gaussian process よりも重い裾を持ちうる sub-Wiebull process に対する Lasso型推定量の誤差評価と，Gaussian process に対する l0 正則化推定量の誤差評価についても述べ，i.i.d. case との比較を行った．さらに，上記の推定量の漸近的な性質を議論する際の設定について詳細な議論を行った．

References
van de Geer, S. A. (2016). Estimation and Testing under Sparsity. Springer.

報告書
赤間陽二

2022 年 11 月 9 日“Unbounded largest eigenvalue of large sample correlation matrices: Asymptotics and applications”(大規模標本相関行列の非有界な最大固有値の漸近論と応用) の題目で講演し, 標本相関行列や標本分散行列の最大固有値の極限での振る舞いについて議論した. ただし, 母集団分布は多次元正規分布だが一定の値 ρ があって異なる変数の組みはどれも相関係数が ρ ≥ 0 になっていることを要求した. この母集団を equi-correlated normal population (ENP) ということにした. ENP は, intraclass covariance matrix [3] を母分散行列とする
正規母集団である. 講演を次のとおり行なった:

1. ENP の標本相関行列と標本分散行列の極限固有値分布は, ρ = 0 の場合の極限固有値分布が 1 − ρ 倍縮小されるという, 講演者たちの結果 [1]を説明した.

2. 講演者たちの結果 [1] の拡張として, 一般の ENP に対して, GuttmanKaiser criterion が削減する次元の割合は, 次元 p と標本の大きさ n が共に大きいときに, さらに, それらの比 p/n が非常に小さいときは, ρ が0 ならば 1/2 に収束し, 同相関係数が正ならば 0 に収束し, 母分散は任意であることを注意した. これは, 自由確率論でとりあえず示せること
より強い結果であった.

3. 講演者の結果として, 標本相関行列 R の最大固有値 λ1 (R) を R のサイ
ズ p で割ったものが ρ の強一致推定量であることを紹介した.
4. 経済学のデータとして S&P500 の株のリターンの平均相関係数の時系
列を, 相関構造が equicorrelation である GJR GARCH [2] で計算し, その時間平均 ρ が λ1 (R)/p よりやや小さいが, 両者の間の重回帰係数が高いため, 高速に計算できる λ1 (R)/p が有用であることを指摘した.

5. 分子生物学における binary multiple sequence alignment のデータが,我々の数学的結果 [1] と講演者の結果に比較的合うことを説明した.

6. DNA マイクロアレイデータの相関行列のヒートマップを提示し, それらさまざまな構造を持つにも拘らず, ガットマン・ケイザー基準が示唆する次元と標本の大きさとの積は, ほぼ 1 であることを指摘した。

1

7. 将来の研究方向として, A. Quadeer (香港科学技術大学), D. MoralesJimenez (Queen’s University Belfast), M. McKay (University of Melbourne) のグループとの議論で彼らが示唆していた相関行列のモデル,
すなわち, 生物学や経済学に現れる相関行列で, ブロック構造とバックグラウンドの同相関があり, ただし, ブロックサイズはそんなに大きくならないけれど, ブロックの個数が増えていくモデルを示した.

ENP の標本分散行列の最大固有値の漸近正規性について, 最大固有値のデータ次元に関するオーダーについて注意を喚起され, 青嶋誠氏と矢田和善氏から関連文献 [4, 5, 6] および母集団分布が ENP であるかの検定手順として文献 [3] の付録を紹介された.
最大固有値の次元に関するオーダーについての矢田・青嶋の一連の論文により, 2022 年の秋の数学会 (北海道大学) での講演者の発表に対する, 確率論研究者からの質問に答える見込みを得た. すなわち, ρ > 0 である正規母集団の標本分散行列の最大固有値の漸近正規性についての着想を得た.
高次元 PCA を観測天文学のデータに適用する天文学者に, 高次元の探索的因子分析の可能性について質疑して, 面識を得た.

参考文献[1] Y. Akama and A. Husnaqilati. A dichotomous behavior of GuttmanKaiser criterion from equi-correlated normal population. J. Indones.
Math. Soc., 2022. To appear. arXiv:2210.12580v2.
[2] R. Engle and B. Kelly. Dynamic equicorrelation. J. Bus. Econ. Stat.,30(2):212–228, 2012.
[3] Aki Ishii, Kazuyoshi Yata, and Makoto Aoshima. Hypothesis tests forhigh-dimensional covariance structures. Annals of Institute of Statistical
Mathematics, 2021.
[4] Kazuyoshi Yata and Makoto Aoshima.

PCA consistency for non-gaussian data in high dimension, low sample size context. Communications in Statistics - Theory and Methods, 38(16), 2009.
[5] Kazuyoshi Yata and Makoto Aoshima.

Eﬀective PCA for high-dimension, low-sample-size data with singular value decomposition ofcross data matrix. Journal of Multivariate Analysis, 101(9):2060–2077,
2010.
[6] Kazuyoshi Yata and Makoto Aoshima. Eﬀective PCA consistency forhigh-dimension, low-sample-size data with noise reduction via geometricrepresentations. Journal of Multivariate Analysis, 2012.

2

GMANOVA モデルにおける新たな推定方法と解釈 (報告書)中京大学 教養教育研究院 永井 勇
本講演では, n 個の各個体に対して, 全ての個体で測定時点を揃えて p 回測定して得られる
経時測定データの分析を考えた．このようなデータはバランス型経時測定データと呼ばれ, 各個体で測定時点が揃っていないものはアンバランス型経時測定データと呼ばれる. これらのデータの分析の目的は, データの裏に潜む経時変動を上手く捉えることある. 本講演では, バ
ランス型経時測定データの分析について考えた.
バランス型経時測定データの分析の際には, Pothoﬀ and Roy (1964) で提案された次の一
般化多変量分散分析 (Generalized Multivariate Analysis of Variance; GMANOVA) モデルがよく使われるため, 本講演でもこのモデルを用いた;

Y = 1n µ′ X ′ + AΞX ′ + E,

(1)

ここで, 1n は n 次元の全てが 1 からなるベクトル, 0r は r 次元の全てが 0 からなるベクトル,

Y は各行が各個体で測定して得られる経時測定データからなる n × p 行列, A は各個体の特徴を表す測定時点に無関係な k 個の変数からなる rank(A) = k の n × k 行列とし, A′ 1n = 0k

(各説明変数で中心化されている) を満たしているとし, X は後述のように各行が測定時点の関数からなる p × q 行列であり, これらの Y , A, X は既知である. また, µ は q 次元未知ベク
トル, Ξ は k × q 未知行列であり, E は E[E] = 0n 0′p , Cov[vec(E)] = Σ ⊗ In の n × p 誤差行列とし, Σ は正則な p × p 未知行列とする. このモデルにおいて, E[Y ] = 1n µ′ X ′ + AΞX ′の部分が, 経時測定データの分析で目的としている経時変動に対応していることを報告した.
q−1

ここで測定時点を t1 < t2 < · · · < tp として, X の i 列目を (t0i , t1i , . . . , ti

) とすることは, 経時変動を測定時点の (q − 1) 次多項式で推定することに対応することを報告した.
このモデル (1) において, 未知の µ と Ξ の推定としてよく使われる推定量は, 次のリスクを最小にする µ と Ξ を求めることで得られることを報告した;
{(
(
)′ }
)
.
R(µ, Ξ|Σ) = tr Y − 1n µ′ X ′ − AΞX ′ Σ−1 Y − 1n µ′ X ′ − AΞX ′(2)

実際にこのリスクを最小にする µ̂Σ と Ξ̂Σ を求めると, (X ′ Σ−1 X)µ̂Σ = X ′ Σ−1 Y ′ 1n /n,Ξ̂Σ X ′ Σ−1 X = (A′ A)−1 A′ Y Σ−1 X の解となる. これらの方程式を解き µ̂Σ や Ξ̂Σ を得るた
めに, モデル (1) において rank(X) = q を仮定することが多い. この仮定は経時変動の推定に用いる関数を制限しているとも考えられることを報告した.
もし rank(X) = q ならば (X ′ Σ−1 X)−1 が存在するので, µ̂Σ = X ′ Σ−1 Y ′ 1n (X ′ Σ−1 X)−1 /nであり, Ξ̂Σ = (A′ A)−1 A′ Y Σ−1 X(X ′ Σ−1 X)−1 となる. 実際に経時変動を推定する際は Σが未知のため, その不偏推定量 S = Y ′ {In − 1n 1′n /n − A(A′ A)−1 A′ }Y /(n − k − 1) が代わりに使われる. このとき, E[µ̂S ] = µ, E[Ξ̂S ] = Ξ であることを報告した.
一方で, rank(X) < q の場合, (X ′ Σ−1 X)−1 が存在しないため, これらの推定量は得られな
い. そこで本講演では, 永井 (2021) のアイデアと同様に, q1 + · · · + qr = q となる正の整数をqi とし, µ = (µ1 , . . . , µr )′ (µi ; qi 次元ベクトル), X = (X1 , . . . , Xr ) (Xi ; rank(Xi ) = qi の
p × qi 行列), Ξ = (Ξ1 , . . . , Ξr ) (Ξi ; k × qi 行列) として, モデル (1) を次のように書き換えた;

Y =

r
∑

1n µ′i Xi′ +

i=1

r
∑

AΞi Xi′ + E.

(3)

i=1

ここで Xi の各列が 0p でない際に, r = q (つまり 1 列ずつ分ける) とすると, rank(Xi ) =qi = 1 となる. つまり, Xi の各列が 0p でない場合まで考えることができることを報告した.
このように書き換えたモデルにおいても, リスク (2) と同様のリスクを最小にするようなµ1 , . . . , µr や Ξ1 , . . . , Ξr を求めて並べれば, rank(X) < q でも µ および Ξ の推定量が得られるというのが, 本講演のアイデアであった. 実際, リスク (2) の µ と Ξ と X に分割した形のものを代入すれば, モデル (3) に対応したリスク R′ が得られる. R′ を最小にする µ̂i·Σ と Ξ̂i·Σ

(i = 1, . . . , r) として, それぞれの推定量がどのように求まるかを報告した. ここで, これらを求めて (µ̂1·S , . . . , µ̂r·S ) や (Ξ̂1·S , . . . , Ξ̂r·S ) とすると, µ と Ξ の推定量が得られ経時変動の推定ができる. しかしながら, そのまま展開するとそれぞれ求まりそうにないため, 本講演では
永井 (2021) と同様に, A′ 1n = 0k に注意しつつ, 展開を工夫した.
その結果, R′ = (µ1 , . . . , µr に関連する項) + (Ξ1 , . . . , Ξr に関連する項) と分割できることを報告した. したがって, 第一項を最小にする µ1 , . . . , µr を求め, 第二項を最小にする Ξ1 , . . . , Ξrを別々に求めればよい. 実際, 各 ℓ でそれぞれ求めると以下のようになることを報告した;
µ̂ℓ·Σ = (Xℓ′ Σ−1 Xℓ )−1 Xℓ′ Σ−1 Ξ̂ℓ·Σ = (A′ A)−1 A′ Y −
ここで,

∑
j<1

Xj µ̂j·Σ = 0p ,∑∑Y ′ 1n ∑−Xj µ̂j·Σ  ,nj<ℓΞ̂j·Σ Xj′  Σ−1 Xℓ (Xℓ′ Σ−1 Xℓ )−1 ,j<ℓΞ̂j·Σ Xj′ = 0k 0′p である. これは, ℓ = 1 の場合の両方が陽に求j<1

まり, それを代入することで ℓ = 2 のときの推定量が得られることを示していることを報告した. これらにおいて Σ の部分を S に置き換えることで推定量が構築できる. この推定手法は
永井 (2021) と同様に, 何らかの罰則などを用いていないため最適化のための反復計算が不要であることを報告した. さらに, rank(X) に関する制約なしで推定が可能となるので, 経時変動の推定に用いる関数への制約がほぼ全てなくなることを表していることを報告した. また,これらの推定量を用いるとモデル (3) へ入れることで, 経時変動が推定できる.
これらの推定量の解釈などについて, 当日の講演で報告した. また, この推定量の問題点や
永井 (2022) と同様に別の形の推定量についても当日の講演で報告した.
引用文献:[1] Pothoﬀ, R. F. & Roy, S. N. (1964) A generalized multivariate analysis of variance model

useful especially for growth curve problems. Biometrika, 51, 313–326.
[2] Nagai, I. (2011) Modiﬁed Cp criterion for optimizing ridge and smooth parametersin the MGR Estimator for the nonparametric GMANOVA model. Open Journal of
Statistics, 1, 1–14.
[3] 永井 勇 (2021) 高次元小標本における多変量線形回帰モデルでの推定法, 2021 年度統計関連学会連合大会[4] 永井 勇 (2022) 説明変数がランク落ちしている状況での多変量線形回帰における不偏推
定量, 2022 年度統計関連学会連合大会

ロバストダイバージェンスを用いたベイズ推論について
中川 智之 1
1

東京理科大学理工学部情報科学科

ベイズ統計において, 通常の事後分布は外れ値に対しては頑健ではないため, 外れ値の問題は古くから研究されてきた. 近年では正則変動よりさらに裾の重い対数正則変動 (log-Regularly varying)
するモデルに関するベイズ推定も提案され, 頑健性の 1 つである外れ値が含まれたデータを用いた
事後分布と外れ値を除いたデータを用いた事後分布が漸近的に一致すること (Posterior robustness)が示されている (Desgagné, 2015, など).
一方で, 裾の重い分布の構成は位置尺度母数モデルよりも一般的なモデルへの拡張は容易ではなく, それぞれのモデルやデータ構造によって構成が必要になる. そこで, 本研究では, 一般のモデルに拡張が容易であるロバストなダイバージェンスを用いたベイズ推論を考える. ロバストなダイバージェンスを用いたベイズ推定は, Hooker and Vidyashankar (2014) や Ghosh and Basu
(2016), Nakagawa and Hashimoto (2020) などによって提案されている.
本講演では, 以下のダイバージェンス基準 D(y, θ) を用いた一般化事後分布について紹介した.
fθ (y) を仮定した統計モデルとし, y1 , . . . , yn を g(y) からの標本とする. このとき, 一般化事後分
布は,∑π(θ) exp( ni=1 D(yi ; θ))∑nπ(θ | Y ) := ∫Θ π(θ) exp( i=1 D(yi ; θ))dθと定義される. 以下は, それぞれのダイバージェンスに基づく D(y, θ) である.

• Kullback-Leibler divergence: DKL (y; θ) = log fθ (y).
• Density power divergence (Basu et al. (1998)):11
Dα (y; θ) = fθ (y)α −αα+1∫fθ (x)α+1 dx (α > 0).

• γ-divergence (Jones et al., 2001; Fujisawa and Eguchi, 2008):
1
Dγ (y; θ) = fθ (y)γγ{∫fθ (x)

γ+1}−γ/(γ+1)dx−

1(γ > 0).
γ

Density power divergence, γ-divergence に関してはそれぞれロバストのダイバージェンスとして知られており, 様々な観点から外れ値に対する頑健性が示されている. 一方で, Posterior robustnessとの関係は示されていない. 本研究では, Posterior robustness を満たすためのダイバージェンス
基準の条件として以下を与えた.
exp{D(z; θ)}= 1 (∀θ ∈ Θ)
h(z)
|z|→∞
lim

この条件は, lim|z|→∞ fθ (z) = 0(∀θ ∈ Θ) を満たすとき, γ-divergence を用いた一般化事後分布は常に成り立つことが簡単にわかる. 一方で, 通常の事後分布や Density power divergence を用いた
一般化事後分布に関しては, 十分条件を満たさないことが分かった. 図 1 と図 2 はそれぞれの事後
分布のヒストグラムを描いたものである. 図 1 を見てわかる通り, 通常の事後分布は外れ値を含んだ事後分布と外れ値を除いた事後分布は乖離しており, 外れ値に影響を受けることがわかる. 一方で, 図 2 を見ると γ-divergence を用いた一般化事後分布は外れ値を含んだ事後分布と外れ値を除いた事後分布がほぼ同じである.
Localtion parameter

Localtion parameter
3

5

4

2

3

with_outliers
without_outliers

variabledensitydensityvariablewith_outliers
without_outliers

2
1

1

0

0
0.0

0.5

1.0

1.5

2.0

value

図 1: 通常の事後分布

−0.4

0.0

0.4

value

図 2: γ-divergence を用いた一般化事後分布

講演では応用例として, 正規分布の場合と順序付き回帰の場合についての数値実験の詳細を紹介した. 順序付き回帰については, さまざまなリンク関数を柔軟に用いることができ, ロバストなベ
イズ推定を可能にするためには, ロバストなダイバージェンスを用いる必要があることを報告した.

参考文献
Basu, A., Harris, I. R., Hjort, N. L., and Jones, M. (1998). Robust and eﬃcient estimation byminimising a density power divergence. Biometrika, 85(3):549–559.
Desgagné, A. (2015). Robustness to outliers in location–scale parameter model using logregularly varying distributions. The Annals of Statistics, 43(4):1568–1595.
Fujisawa, H. and Eguchi, S. (2008). Robust parameter estimation with a small bias against
heavy contamination. Journal of Multivariate Analysis, 99(9):2053–2081.
Ghosh, A. and Basu, A. (2016). Robust bayes estimation using the density power divergence.
Annals of the Institute of Statistical Mathematics, 68(2):413–437.
Hooker, G. and Vidyashankar, A. N. (2014). Bayesian model robustness via disparities. Test,
23(3):556–584.
Jones, M., Hjort, N. L., Harris, I. R., and Basu, A. (2001). A comparison of related density-basedminimum divergence estimators. Biometrika, 88(3):865–873.
Nakagawa, T. and Hashimoto, S. (2020). Robust bayesian inference via γ-divergence. Communications in Statistics-Theory and Methods, 49(2):343–360.

代表点を用いた大規模クラスタリングの近似法とその性質
寺田吉壱 1,3 , 山本 倫生 2,3
1 大阪大学大学院基礎工学研究科, 2 大阪大学大学院人間科学研究科
3 理化学研究所革新知能統合研究センター

1. はじめに
近年，データの大規模化と複雑化が進み, データから仮説や有益な情報を獲得することが課題となっており, 探索的なデータ解析の中でも教師なし学習の重要性が再認識されている．クラスタ
リング法は, データの背後のクラスタ構造を明らかにするための教師なし学習の方法であり, 様々な分野で広く応用されている. 最も代表的なクラスタリング法として, k-means 法が挙げられる.
k-means 法はその簡便性と計算コストの低さから多用されるが, その単純さ故にデータの背後にある複雑なクラスタ構造を十分に捉えられない可能性がある. そのため, spectral clustering (vonLuxburg, 2007) など, より柔軟にクラスタ構造を捉えられる方法の利用が望ましい. しかし, これらの複雑な方法は, 一般に計算コストが高く, 計算コストの削減が大きな課題となっている.
カーネル法におけるクラスタリング法に対しては, カーネル法に特化した Nyström 近似や Random
Fourier Feature などの計算コスト緩和法が適用出来る. 一方で, Yan et al. (2009) では, spectralclustering の k-means 法に基づく近似法 (KASP) が提案されている. KASP では, クラスタ数を多く設定した k-means 法を大規模データに適用し, 得られたクラスタ中心をデータを代表する点とする. そして, 得られた代表点のみに spectral clustering など複雑なクラスタリング法を適用し,代表点に対するラベルをその代表点に近いデータ点のラベルとする. KASP は, spectral clusteringに限らず任意のクラスタリング法に対して適用することができる. この方法の大きな利点は, 安定性, 簡便性, 計算コストの低さ, 汎用性である. また, KASP のアルゴリズムとは本質的に異なるが, 代表点として subsample を用いる spectral clustering の近似法が提案されている (Mohan andMonteleoni, 2017). しかし, Terada and Yamamoto (2019) の理論から, この方法は subsample に対する normalized cut と等価であり, 汎用的な近似方法ではない.
本発表では, KASP の問題点を明らかにし, その問題点を解決した汎用的な大規模クラスタリン
グの近似法を提案する.
2. KASP の問題点と提案手法
KASP では, データの代表点として, k-means 法のクラスタ中心を用いる. しかし, k-means 法によって生成した代表点から構成される経験分布は, 母集団分布を代表する点とはならない. 具体
的には, 代表点の経験分布は, 母集団分布よりも裾が重い分布に収束することが示せる. そのため,
KASP を用いて近似を行うと, グラスタリング結果にズレが生じてしまう.
この問題点の最もシンプルな解決法は, 各代表点に適切な重みを与えることである. K-means 法に対応するベクトル量子化は, 母集団分布との L2 -Wasserstein 距離を最小にするような #(supp(Q)) ≤K を満たす離散測度を求める問題に対応している. このことから, 代表点に適切な重みを与えることで, KASP の問題点を解消することができる. 一方で, 代表点の経験分布のズレから, 母集団分
布において密度の低い点も代表点として生成される. それらの代表点には, 低い重みが割り振られるため, 効率が悪い.
この問題を解決するために, 新しい代表点の生成方法である Density-Preserving Vector Quantization (DPVQ) を提案する. DPVQ は重み付き k-means 法の一種であり, 容易に代表点を生成できる. また, DPVQ が生成する代表点の経験分布は, 漸近的にデータの背後の分布へ収束すること

Algorithm 1 VQn (µ | r, K) の最適化アルゴリズム(0)

1: t ← 0 とし, クラスタ中心 µ1(0), . . . , µK を初期化する．

2: for t = 0, . . . , T do
3:(t)

各 i (i = 1, . . . , n) に対して， xi − µk()(t)
U (t) = uijを得る.

を最小にするクラスタ k を割り当て, 帰属行列n×K
1 if ∀j; x − µ(t) ≤ x − µ(t) ,
i
i
(t)
j
k
uik =
0 otherwise.
4:

クラスタ平均を以下で更新する.
(t+1)

µ̂k

=∑

1

n
∑

(t) (t)n
j=1 ujk wjk i=1

(t)

(t)

uik wik xi ,

(t)

wik =


 x − µ̂(t) r−2
i

δ

k

if xi − µ̂k

(t)

> 0,

(t)
xi − µ̂k= 0.

if

ここで, δ > 0 は小さい正の定数である.
5:
収束判定条件を満たせば停止し, 満たさなければ t ← t + 1 とする.
6: end forが示せる. そのため, DPVQ による代表点には平等な重みが割り振られるため, 効率的な近似が期待できる. 一方で, DPVQ は, 密度推定を必要とするため, 高次元データに対しては不安定となる.
そこで, 本発表では, 以下で定義される order r のベクトル量子化器を用いた近似法も提案する.

1∑
min ∥xi − µk ∥r
VQn (µ | r, K) :=
1≤k≤K
n
n

i=1

ここで, r ∈ (0, 2] は定数, ∥ · ∥ は Rd 上のノルム, x1 , . . . , xn ∈ Rd は各データ点, µk ∈ Rd は k 番目のクラスタ中心, µ = (µ1 , . . . , µK ) である. 本発表では, 計算の簡便性のために, ∥ · ∥ は Euclid ノ
ルムとする. r = 2 とすれば VQn (µ | r, K) は k-means 法と一致するが, r を小さくとることで代表点の経験分布と母集団分布のズレを小さくすることができる. 本発表では, VQn (µ | r, K) に対する最適化問題を高速に解くために, k-means like なアルゴリズムを提案する. Algorithm 1 において, VQn (µ(t) | r, K) > VQn (µ(t+1) | r, K) という単調減少性が成り立つため, 停留点への収束性が保証できる. 提案手法を用いた場合の spectral clustering の近似方法, 詳細な理論的性質, DPVQ
や r < 2 としたベクトル量子化器を用いた提案手法と KASP や subsample を用いた既存の近似手
法の数値実験による比較は当日報告する.
参考文献

Mohan, M. and Monteleoni, C. (2017). Beyond the Nystrom Approximation: Speeding up
Spectral Clustering using Uniform Sampling and Weighted Kernel k-means. In Proceedingsof the 26th International Joint Conference on Artificial Intelligence, 2494–2500.
Terada, Y. and Yamamoto, M. (2019). Kernel Normalized Cut: a Theoretical Revisit. In
Proceedings of the 36th International Conference on Machine Learning, PMLR 97, 6206–
6214.
von Luxburg, U. (2007). A tutorial on spectral clustering. Statistics and Computing, 17, 395–
416.
Yan, D., Huang, L., and Jordan, M. I. (2009) Fast approximate spectral clustering. In Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data
Mining, 907–916.

Recent progress in the Application of High-Dimensional Statistics to
Astrophysics and Cosmology
竹内 努 1,2, COORAY, Suchetha1, 矢田 和善 3, 青嶋 誠 3, 石井 晶 4, 江頭 健斗 5,
吉川耕司 6, 中西 康一郎 7, 河野 海 1, 河野 孝太郎 8
1. 名古屋大学理学研究科素粒子宇宙物理学専攻,
2. 統計数理研究所統計的機械学習研究センター,
3. 筑波大学数理物質系,
4. 東京理科大学理工学部情報科学科,
5. 筑波大学 理工情報生命学術院,
6. 筑波大学計算科学研究センター,
7. 国立天文台 ALMA プロジェクト,
8. 東京大学天文学教育研究センター
銀河とは, 星と星間物質(ガスと星間塵), 暗黒物質からなる巨大な天体である. 銀河は現在観測可能な宇宙に数千億個存在しており, 我々の目に見える波長(可視光線)での宇宙の姿を形作っている. しかし宇宙誕生当時, 物質分布はほぼ一様であった. 銀河は平均からわずかに密度が高い領域が重力によって成長し, 合体成長を経て現在の姿へと進化してきたのである. 銀河進
化は周囲の銀河の密度やガス密度など, 銀河の置かれた環境にも強く依存する極めて複雑な過程である. そして銀河の形成と進化は 138 億年の宇宙進化の歴史の中でも非常に重要な現象のひとつであると考えられている.
実験室実験の不可能な遠方宇宙を対象とする天文学において, 天体の分子の物理情報を得る唯一の手段が分光観測である. しかし詳細な分光観測は一般的に非常に時間がかかり, 天体をマッピングして独立な観測点を多数得ることは容易ではない. また, 対照が非常に稀な天体であることも多い. どちらのケースも, データの持つ波長方向の次元 d に対してサンプル数 n が非常に小さい(n ≪ d)という状況が生じる. 伝統的な天文学では, このような状態は不良設定問題とされ, 物理的考察のためにはデータ d の情報を大幅に捨て, n より小さくして解析する他ないと考えられていた. しかし当然, このような無駄は避け, 情報量をフルに活かせる方法が望ましい. n ≪d となるデータは高次元小標本(HDLSS)データと呼ばれる. 天文学以外の分野, たとえばゲノム
解析では遺伝子の標本数 n ~ 100 に対し塩基配列の次元数が d ~ 105 といった問題は特殊ではない. HDLSS の問題を解決するためにここ 10 年で発展した統計学の方法が高次元統計学であり, 現在も続々と新しい知見が得られている(e.g., Aoshima 2018).
本研究では, 1) 電波分光マップ観測, および 2) 希少天体の分光観測の 2 つの場合について,
高次元統計の方法を適用し, 従来の天文学では不可能とされてきた解析を行う. 1)では, 爆発的星形成銀河 NGC253 の各点をチリの大型ミリ波サブミリ波干渉計 ALMA によって分光観測した探査(マッピング)型について解析する. 分光マッピングには一般的に大変な観測的コストがかかり,
ALMA の場合現状では最大級でも n ~ 200 程度の観測点が得られる程度であるのに対して波長(振動数)分散方向には d ~ 2000 の情報がある. よって, ALMA の分光マッピングは典型的な

HDLSS データである. 従来の天文学では, 分光データから物理の第一原理から有用であるとわかっているいくつかの輝線のみを抽出し, それを用いて分類を行っていた. しかし, ALMA の観測によって, たとえ隣り合った星間分子雲の領域であっても分子輝線は大幅に異なった様相を示すことが判明し, 古典的方法の限界が明らかになっている. 巨大すぎる情報量のため, 分子雲の進化と星形成の関係が全く見えないのである. 本研究では NGC253 の ALMA マップに高次元主成
分解析(PCA)を適用した. この結果, 高次元 PCA は物理的仮定をおかずに NGC253 中心部のガスの回転運動や星形成領域からのガス流(アウトフロー)を自動的に検出できることを発見した.
これは, 複雑に見える分光データはわずか数個の固有スペクトルによって表現できることを意味
する. Yata & Aoshima (2022)の方法によって, それぞれの固有スペクトルをコントロールするのはHCN (青酸)および HNC 分子であることも見出した. 今回の解析は比較的振動数分解能の小さい
データについて行ったが, よりよいデータがすでに得られており, 星形成の物理に肉薄できる期待が持てる.
さらに 2)のタイプとして, クェーサーの原子水素吸収線系(HI forest)による「宇宙再電離期」の検証を進めている. 宇宙誕生から 1 億年以内の初期, 宇宙は原子状態の水素で満たされていた状態から電離水素がほとんどの体積を占める状態への相転移を経験した. これが宇宙再電離と呼ばれている現象である. この研究に有効なのが, 宇宙初期に形成されたクェーサーと呼ばれる天体である. クェーサーは巨大なブラックホールに物質が落下することにより, 重力エネルギーが解放されて巨大なエネルギーの放射が生じている天体で, あらゆる電磁波の波長できわめて明るい.
我々は長波長電波(波長数 cm－数 m)の放射に注目した. クェーサーの手前に水素ガスが存在すると, クェーサーの明るい電波放射に吸収線という「影」として見える. この吸収線の数は膨大であり, 天体形成および宇宙論について豊富な情報を持っていることが知られている. しかし, 宇宙
再電離期にすでに形成され, 活動していたクェーサーの数はかなり少ないことが明らかになっており, 現在建設中の次世代超大型長波長電波干渉計 Square Kilometre Array (SKA)をもってしても数十個しか検出されないと予想されている. すなわち, クェーサーの HI forest の観測も典型的な HDLSS データである. 我々は HI forest から空間 2 点相関関数およびその分散, そして宇宙
年齢依存性を高次元統計の方法で解析する方法を開発した. 宇宙論シミュレーション Illustris
TNG のデータを用いた検証から, この新しい手法によって HI forest を通じた宇宙再電離期の物
質の空間分布および進化を検出できることが見出された. 今後の発展が期待される.

凸クラスタリングによるマルチタスク学習
電気通信大学 大学院情報理工学研究科
九州大学 大学院数理学研究院

岡崎
川野

彰良
秀一

1 はじめにマルチタスク学習とは，タスクと呼ばれる複数のデータ集合が得られている際にタスク毎にモデルを設定し，タスク間の関係性を考慮しながらモデルの推定を行う方法論である．タスク毎に設定された各モデル
が他のモデルの情報を取り入れることで，モデル全体の推定精度が改善される．ネットワーク lasso 正則化(Hallac et al., 2015) を用いた手法 (Yamada et al. (2017); He et al. (2019)) では，モデル毎に設定された回
帰係数ベクトルの値を縮小して推定することにより，クラスタリングに基づくマルチタスク学習を実行している．しかし，これらの手法の問題点として，ネットワーク lasso 正則化は推定量としてバイアスを持ち，異なる
クラスタに属するタスクの情報に影響を受けやすい問題が存在する．そこで本報告では，ネットワーク lasso
正則化に基づくマルチタスク学習の枠組みにおいて，クラスタの重心を表す新たなパラメータを各タスクに対して導入したマルチタスク学習手法を提案する．クラスタの重心に関するパラメータを導入したことにより，縮小によってバイアスを受ける推定量と回帰を行うための推定量が分離され，推定のバイアスが軽減される．
提案手法は k 平均法に基づくマルチタスク学習手法の Zhou et al. (2011) とは異なる形での凸緩和手法として見なすことができ，解釈が容易である．
提案手法の有効性をモンテカルロ・シミュレーションおよび実データに対する適用を通して検証する．

2 グループ連結正則化に基づくマルチタスク学習
nm 個の説明変数に関する p 次元データ {xmi ; i = 1, . . . , nm } および nm 個の目的変数に関するデータ{ymi ; i = 1, . . . , nm } に対して，これらの組 {(ymi ; xmi ), i = 1, . . . , nm } が独立に与えられているとし，目的
変数ベクトル ym = (ym1 , . . . , ymnm )⊤ ∈ Rnm および計画行列 Xm = (xm1 , . . . , xmnm )⊤ ∈ Rnm ×p として表す．これらの組を一つのタスクとし，さらに T 個の関連するタスク {(ym , Xm ); m = 1, . . . , T } が与えられているとする．T 個のタスクに対して，目的変数と説明変数の関係を次の T 個の線形回帰モデル

ym = Xm w m + ϵm ,m = 1, . . . , T(1)により捉える．ここで，wm = (wm1 , . . . , wmp )⊤ は m 番目のタスクに対応する回帰係数ベクトルであり，ϵmは各要素がそれぞれ独立に N (0, σ 2 ) に従う観測誤差ベクトルである．(1) 式に対して，次の最小化問題


T
∑



∑
1
2
min
∥y
−
X
w
∥
+
λ
r
∥w
−
w
∥
m
m m 2
m,l
m
l q
wm ∈Rp , m=1,...,T 

2nm
m=1

(2)

(m,l)∈Eを考える．ここで，rm,l は非負の値を取る m 番目と l 番目のタスクの関係の大きさを表す重みであり，E はタ
スクの組の集合を表す．また，λ は非負の値を取る正則化パラメータである． 第二項目は m 番目と l 番目のタスクに関する回帰係数ベクトル wm , wl の併合を誘引する Lq ノルムグループ連結正則化項である．(2) 式は1

q ≥ 1 である場合，凸最適化問題であることが知られており，大域的最適解が得られる．また，q = 2 のとき，
ネットワーク lasso 正則化 (Hallac et al., 2015) の枠組みとして捉えることができる．

3 提案手法: 凸クラスタリングに基づくマルチタスク学習
最小化問題 (2) においてクラスタリングを実行する第二項目の正則化項は，異なるクラスタに属する回帰係数ベクトル間で縮小を行い，推定精度が低下する問題を含んでいる．この問題を軽減するため，Yamada et al.

(2017) や He et al. (2019) では凸クラスタリングと同様に k-近傍法によって重み rm,l の値を与えている．また， Zhou and Zhao (2016) や Shimamura and Kawano (2021) では rm,l をモデルに含まれる潜在変数として考え，その計算アルゴリズムを提案している．しかし，最小化問題 (2) における重み rm,l を潜在変数として推定する方法は目的関数の非凸性を誘引し，大域的最適解を得ることが困難になる問題が生じる．
この回帰係数ベクトル間の縮小に関する問題を克服するために，本報告では以下の最小化問題を考える．


T
∑


T

∑
∑
1
λ
1
2
2
min
∥y
−
X
w
∥
+
∥w
−
u
∥
+
λ
r
∥u
−
u
∥
.
m
m
m
m
m
2
m,l
m
l
2
2
2
wm ,um ∈Rp , m=1,...,T 

2nm
2 m=1
m=1
(m,l)∈E(3)
ここで，um ∈ Rp は m 番目のタスクが属するクラスタの重心に関するパラメータ，λ1 , λ2 は正の値を取る正則化パラメータである．本報告では，この手法を MTLCVX (Multi-Task Learning via ConVeX clustering)と呼ぶ．
最小化問題 (3) の第二項目および第三項目は wm に関して凸クラスタリングを実行する項である．最小化問題 (2) とは異なり，回帰係数 wm に関しては直接縮小を行わず，um に関して縮小を行う．そして，第二項目により，縮小推定された um の周りで wm を推定することでマルチタスク学習が実行される．これにより，

wm に関して生じる縮小に関するバイアスが軽減される．また，最小化問題 (3) は um , wm に関して同時凸関数であるため，その解は大域的最適解となる．モンテカルロ・シミュレーションおよび実データ解析の結果については当日報告する。

参考文献
Hallac, D., Leskovec, J., and Boyd, S. (2015). Network lasso: Clustering and optimization in large graphs.
In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data
Mining, 387–396.
He, X., Alesiani, F., and Shaker, A. (2019). Eﬃcient and scalable multi-task regression on massive
number of tasks. In Proceedings of the AAAI Conference on Artiﬁcial Intelligence, 33(01), 3763–3770.
Shimamura, K. and Kawano, S. (2021). A bayesian approach to multi-task learning with network lasso.
Preprint, arXiv:1402.6455.
Yamada, M., Koh, T., Iwata, T., Shawe-Taylor, J., and Kaski, S. (2017). Localized Lasso for HighDimensional Regression. In Proceedings of the 20th International Conference on Artiﬁcial Intelligenceand Statistics, 54, 325–333.
Zhou, J., Chen, J., and Ye, J. (2011). Clustered multi-task learning via alternating structure optimization.
Advances in Neural Information Processing Systems, 24, 702–710.
Zhou, Q. and Zhao, Q. (2016). Flexible clustered multi-task learning by learning representative tasks.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(2), 266–278.

2

馬蹄事前分布に基づく連結 lasso ロジスティック回帰モデリング1

電気通信大学大学院情報理工学研究科

書川侑子

九州大学大学院数理学研究院

川野

秀一

はじめに

二値データの解析手法として，ロジスティック回帰モデルに正則化項を取り入れること
でデータのスパースな構造を予測に反映させる手法は近年注目を集めている．ベイズ統
計学の観点に立つと，正則化は回帰係数に縮小事前分布を仮定することとして捉え直すことができる．中でもベイジアン連結 lasso (Kyung et al., 2010) は，回帰係数と隣接する回帰係数の差にラプラス事前分布を仮定することで，変数選択および変数の併合を可能にする．しかし，ラプラス事前分布には本来非ゼロと推定されるべき箇所を過剰にゼロに縮小するという問題点がある．
この問題点を解決するため，本報告では，馬蹄事前分布 (Carvalho et al., 2010) を連続する回帰係数の差に仮定したロジスティック回帰モデルを提案した．推定方法としては
Polson et al. (2013) による尤度の階層表現を用いることでギブスサンプリングを可能にする．また，数値実験を通して提案手法の有効性を検証した．

2

馬蹄事前分布による連結 lasso を用いたロジスティック回帰モデル

二値を取る目的変数 y1 , y2 , · · · , yn と p 次元説明変数ベクトル x1 , x2 , · · · , xn に対し, 次の二値ロジスティック回帰モデルを考える．

Pr(yi = 1 | xi ) =1
1 + e−(β0 +xi β)T,

i = 1, . . . , n.

(1)

ここで，Pr(yi = 1 | xi ) は xi が与えられたもとで yi = 1 となる確率，β0 は切片項，
β = (β1 , β2 , · · · , βp )T は回帰係数ベクトルである．
(1) 式のロジスティック回帰モデルにおいて，回帰係数にラプラス分布を仮定し，回帰係数の差に馬蹄事前分布を仮定することを考える．Pólya-Gamma 分布を用いた尤度の階層表現 (Polson et al., 2013) を用いて，提案モデルは次のように表される．
()
1yi | xi , β, β0 , wi ∼ Binom 1,,T
1 + e−(β0 +xi β)

wi ∼ PG(1, 0),()βj | λ̃1 ∼ Laplace 0, λ̃1 ,()βj − βj−1 | λ2j , τ̃ 2 ∼ N 0, λ2j τ̃ 2 ,λj ∼ C+ (0, 1),τ̃ ∼ C+ (0, 1).

1(2)

ここで，パラメータ a > 0, b を持つ Pólya-Gamma 分布は次の確率密度関数で表される．
∞
1 ∑gk
PG(x | a, b) = 21c2
2π k=1 (k − 2 )2 + (4π2)

gk はガンマ分布 Ga(a, 1) に従う独立な確率変数である．また，Laplace(x | 0, λ) は，確率λ密度関数 Laplace(x | 0, λ) = exp (−λ|x|) で表されるラプラス分布であり，C+ (x | 0, 1)22+は確率密度関数 C (x | 0, 1) =で表される半コーシー分布である. 馬蹄事前分
2
π (x + 1)
布を連続する回帰係数の差に仮定することで，回帰係数の差のうち非 0 に推定されるべき箇所が過度に縮小されることを抑えた推定が行える．
正規尺度混合 (Andrews and Mallows, 1974) と半コーシー分布の階層表現 (Wand et al.,2011) を用いると (2) 式のモデルは次のように表せる.
()
1
yi | xi , β, β0 , wi ∼ Binom 1,,T
1 + e−(β0 +xi β)
wi ∼ PG(1, 0),βj | τj2 ∼ N(0, τj2 ),( )λ̃21τj2 ∼ EXP
,
2
(
)
1 1
2
,
λj | νj ∼ IG
,
2 νj
)
(
1 1
2
,
,
τ | ξ ∼ IG
2 ξ
(
)
1
ν1 , . . . , νp , ξ ∼ IG
,1 .
2

(3)

(3) 式の階層表現により，回帰係数の推定方法としてギブス・サンプリングが構成できる．

参考文献
Andrews, D. F. and Mallows, C. L. (1974). Scale mixtures of normal distributions. Journalof the Royal Statistical Society: Series B (Methodological), 36(1), 99–102.
Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010). The horseshoe estimator forsparse signals. Biometrika, 97(2), 465–480.
Kyung, M., Gill, J., Ghosh, M., and Casella, G. (2010). Penalized regression, standard
errors, and Bayesian lassos. Bayesian Analysis, 5(2), 369–411.
Polson, N. G., Scott, J. G., and Windle, J. (2013). Bayesian inference for logistic models
using pólya–gamma latent variables. Journal of the American Statistical Association,
108(504), 1339–1349.
Wand, M. P., Ormerod, J. T., Padoan, S. A., and Frühwirth, R. (2011). Mean ﬁeld
variational Bayes for elaborate distributions. Bayesian Analysis, 6(4), 847–900.
2

凸最適化によるエルゴード的拡散過程のオンライン推定
ONLINE ESTIMATION OF ERGODIC DIFFUSION PROCESSESWITH CONVEX OPTIMIZATION
仲北祥悟（東京大学大学院総合文化研究科）
Let us consider the parametric estimation of the following d-dimensional stochastic differential
equation (SDE):dXta,b = b Xta,b dt + a Xta,b dwt , X0 = x ∈ Rd , t ≥ 0.
SDEs describe dynamics with randomness and allow for flexible model structures under mild
conditions. Therefore, they are used to model phenomena in broad disciplines such as finance,
biology, epidemiology, physics, meteorology, and machine learning. In this study, we proposea,b
an online parametric estimation method of b based on discrete observations {Xih}withn i=0,...,nhn > 0.
Online estimation, where the estimator is updated as data are acquired, is a typical andsignificant concern in time series data analysis because it is quite useful for real-time decision
making. However, most studies on the online parametric estimation of SDEs depend on thesetting of continuous observations {Xta,b }t≥0 [3, 1, 2], which is restrictive in real data analysis.
Hence, we aim to propose online estimation methods for SDEs with discrete observations.
We provide uniform risk bounds for the parametric estimation of both diffusion and drift
coeﬀicients of SDEs with discrete observations and model misspecification via online gradientdescent with convex loss functions and their convex approximations. Those bounds give theoretical convergence guarantees of the proposed online estimation method for SDEs with discreteobservations, which are the main contribution of our study. To derive the bounds, we combine
the three theoretical discussions: (i) model-wise non-asymptotic risk bound for the stochastic
mirror descent (SMD) with dependent and biased subgradients; (ii) simultaneous ergodicity
and uniform moment bounds for a class of SDEs; and (iii) the proposal of loss functions for the
online parametric estimation.
Selecting drift estimation as an example in this section, we set the convex and compact
parameter space Θ ⊂ Rp and the triple of measurable functions (bm , M, J) such that bm (x, θ)is the possibly misspecified parametric model, M (x) is a positive semi-definite weight function,
and J (θ) is the regularization term. We set the functionhi
1
ϕ (x, y, θ) := M (x) (y − bm (x, θ))⊗2 + J (θ) .
2
Assume that ϕ (x, y, θ) is convex in θ for all x, y ∈ Rd and has measurable elements in thesubdifferential for all x, y, and θ. {θi ; i = 1, . . . , n + 1} defined by the following online gradientdescent algorithm



hn
1
a,b
a,b
,
θi+1 := ProjΘ θi − √ ∂θ ϕ X(i−1)h
∆
X
,
θ
,
i
i
n h
i
n
a,b
with an arbitrary initial value θ1 ∈ Θ and a sequence of discrete observations {Xih;i =n0, . . . , n}, is then well-defined as a sequence of random variables by choosing measurable suba,ba,bgradients, where ∆i X a,b = Xih− X(i−1)hand hn > 0 is the discretization step. Note that thennlearning rate chosen here does not lead to the best convergence but is simple and approximately
the best in our study. Our contributions (i) and (iii) provide the following risk bound for theThis work was partially supported by JSPS KAKENHI Grant Number JP21K20318 and JST CREST Grant
Number JPMJCR21D2 and JPMJCR2115.
12

S NAKAKITA

estimator θ̄n := n1Pni=1 θi with a fixed (a, b): for some c > 0,hsup Ea,b
f a,b θ̄nx

i

− f a,b (θ) ≤ cθ∈Θ

log nh2np+ hβ/2nnh2n!
,where β ∈ [0, 1] is a parameter controlling the smoothness of b, Ea,b
x is the expectation overR
{Xta,b }t≥0 with X0a,b = x, f a,b (θ) = M (ξ) [(bm (ξ, θ) − b (ξ))⊗2 ]Πa,b (dξ) + J (θ) is the loss
function, and Πa,b is the invariant probability measure of Xta,b . The contribution (ii) yields theexistence of c such that the inequality holds uniformly in S := {(a, b)}, a class of coeﬀicientsof SDEs satisfying the same regularity conditions; hence, the risk bound is uniform in S. Notethat θ̄n estimates the best θ∈ Θ (or the quasi-optimal parameter; see [4]) with bm (, θ) closest
to the true b in the L2 Πa,b -distance. Moreover, if the model bm correctly specifies b, that is,for all (a, b) ∈ S there exists θ such that b = bm (·, θ), then we can obtain the following bound:
!
h
i
2

log
nh
sup Ea,b
f a,b θ̄n ≤ c p n + hβ/2
.
x
n
nh2n
(a,b)∈S
One simple but significant outcome of the above discussion is a non-asymptotic risk guaranteeof the following online gradient descent for linear models such that an arbitrary initial valueθ1 ∈ Θ,


 

1  m  a,b
a,b
a,b
m
θi+1 := ProjΘ θi + √ ∂θ b X(i−1)hn , θi∆i X − hn b X(i−1)hn , θi,iwhere bm (x, θ) is the possibly misspecified parametric model whose components are linear inθ ∈ Θ. Note that it corresponds toPthe case M (x) = Id , J (θ) = 0. As evident, the uniform riskbound for the estimator θ̄n := n1 ni=1 θi over a certain family S of the coeﬀicients a, b holds:
for some c > 0,

Z
 Z


2 a,b
2 a,b
a,b
m
m
sup sup Exb ξ, θ̄n − b (ξ) 2 Π (dξ) − ∥b (ξ, θ) − b (ξ)∥2 Π (dξ)(a,b)∈S θ∈Θ≤clog nh2np+ hβ/2n
nh2n!
.

If we assume that bm correctly specifies b, then
Z


2 a,b
a,b
m
sup Exb ξ, θ̄n − b (ξ) 2 Π (dξ) ≤ c(a,b)∈S

log nh2np+ hβ/2
n
2
nhn

!
.

References[1] Bhudisaksang, T. and Cartea, A. (2021). Online drift estimation for jump-diffusion processes.
Bernoulli, 27(4):2494–2518.
[2] Sharrock, L. and Kantas, N. (2022). Joint online parameter estimation and optimal sensor
placement for the partially observed stochastic advection-diffusion equation. SIAM/ASA
Journal on Uncertainty Quantification, 10(110):55–95.
[3] Surace, S. C. and Pfister, J.-P. (2019). Online maximum-likelihood estimation of the parameters of partially observed diffusion processes. IEEE Transactions on Automatic Control,64(7):2814–2829.
[4] Uchida, M. and Yoshida, N. (2011). Estimation for misspecified ergodic diffusion processes
from discrete observations. ESAIM: Probability and Statistics, 15:270–290.

NON-SPARSE HIGH-DIMENSIONAL ASYMPTOTICS:
THEORY AND PRACTICEMASAAKI IMAIZUMI1
1 The University of Tokyo

1. Introduction
This talk will deal with the recently developed high-dimensional asymptotics. Traditionally, high-dimensional statistics based on sparsity havebeen mainly investigated. In contrast to the literature, in recent years, several high-dimensional analysis without sparsity has been attracting attention,
motivated by the recent success of non-sparse large-scale statistical analysis. For example, there is a high-dimensional analysis using the spectrum ofdata matrices and an analysis of high-dimensional models using the Gaussian comparison theorem. In this presentation, I will overview those theories
and introduce some research results applying them.
This presentation consists of several studies. First, we present an extension of the concentration inequality for empirical averages of highdimensional random matrices (Koltchinskii and Lounici, 2017), and related
over-parameterized linear regression studies (Bartlett et al., 2020). We extend these results to distributions with thick hems. We then present methods
for dealing with high-dimensional models under Gaussianity. These aredone using the Gaussian comparison theorem (Thrampoulidis et al., 2015,2018) and the proud message propagation method (Bayati and Montanari,
2011; Sur and Candès, 2019). We develop M-estimators under generalregularization as well as generalizations of these estimators in the case of
unknown link functions.
Date: December 8, 2022.
1

2

NON-SPARSE HIGH-DIMENSIONAL ASYMPTOTICS: THEORY AND PRACTICE

References
Bartlett, P. L., Long, P. M., Lugosi, G. and Tsigler, A. (2020) Benign
overfitting in linear regression, Proceedings of the National Academy of
Sciences, 117, 30063–30070.
Bayati, M. and Montanari, A. (2011) The lasso risk for gaussian matrices,
IEEE Transactions on Information Theory, 58, 1997–2017.
Koltchinskii, V. and Lounici, K. (2017) Concentration inequalities and moment bounds for sample covariance operators, Bernoulli, 23, 110–133.
Sur, P. and Candès, E. J. (2019) A modern maximum-likelihood theory
for high-dimensional logistic regression, Proceedings of the National
Academy of Sciences, 116, 14516–14525.
Thrampoulidis, C., Abbasi, E. and Hassibi, B. (2018) Precise error analysisof regularized 𝑚-estimators in high dimensions, IEEE Transactions on
Information Theory, 64, 5592–5628.
Thrampoulidis, C., Oymak, S. and Hassibi, B. (2015) Regularized linear
regression: A precise analysis of the estimation error, in Conference on
Learning Theory, PMLR, pp. 1683–1709.

Sparse Bayesian inference on gamma-distributed observations
Yasuyuki Hamura1 , Takahiro Onizuka2 , Shintaro Hashimoto2 , Shonosuke Sugasawa3
1
Graduate School of Economics, Kyoto University
2
Department of Mathematics, Hiroshima University
3
Center for Spatial Information Science, The University of Tokyo

In various statistical applications, we often face a sequence of positive-valued observations
such as machine failure time, store waiting time, survival time under a certain disease, an
income of a certain group, and so on. A common feature of the data is “sparsity” in the sensethat most of the underlying means of observations are concentrated around a certain value(grand mean) while a small part of the means is signiﬁcantly away from the grant mean. To
reﬂect the sparsity structure, a useful Bayesian technique is an idea of “global-local shrinkage”(e.g. Carvalho et al., 2010; Polson and Scott, 2010; Datta and Dunson, 2016) that provides
adaptive and ﬂexible shrinkage estimation of underlying means; when the observations arearound the grand mean, the posterior mean strongly shrinks the observation toward the grandmean, but the observations that are away from the grand mean remain unshrunk.
In this presentation, we consider the following gamma inverse-gamma hierarchical model:δiyi | λi ∼ Ga δi ,λi

,λi | ui ∼ IG(1 + τ ui , βτ ui ),ui ∼ π(·),

i = 1, . . . .n,

where β > 0 and τ > 0 are unknown hyper-parameters for which we assign conjugate gamma
priors, and λi is a parameter of interest. The prior mean of λi is E[λi ] = β, so that β isinterpreted as a grand mean of underlying heterogeneous means. A ﬁxed constant δi is selected
in given context. As considered in Liu and Stephens (2016), if yi and λi are sampling andtrue variances, respectively, the choice is δi = ni /2, where ni is a sample size used to compute
yi . Moreover, if yi is a sample mean based on ni samples generated from an exponential
distribution Exp(1/λi ), it holds that δi = ni , and it reduces the framework of a sequence of
exponential data when ni = 1, considered in Donoho and Jin (2006). Unlike existing shrinkage
priors, our new prior is a shape-scale mixture of inverse-gamma distributions, which has adesirable interpretation of the form of posterior mean and admits ﬂexible shrinkage. In fact,
the posterior mean of λi is given byδi yi + βτ ui
E[λi | yi ] = E
δi + τ u i


yi = yi − E[κi | yi ](yi − β),where κi = τ ui /(δi + τ ui ) ∈ (0, 1) is known as shrinkage factor that determines the amountof shrinkage of yi toward the grand mean β. As desirable properties of κi , E[κi | yi ] should
be close to 1 when yi is close to the grand mean, leading to strong shrinkage toward β, whileE[κi | yi ] should be suﬃciently small for yi having large yi − β to prevent bias caused byover-shrinkage. We also note that the global parameter τ determines the overall shrinkage
eﬀect, whereas the local parameter ui allows κi to vary over diﬀerent observations. As priorsof ui that satisfy such desirable properties, we propose the following scaled beta (SB) andinverse rescaled beta (IRB) priors.
1
ui a−1,
B(a, b) (1 + ui )a+b1{log(1 + 1/ui )}b−1
1,πIRB (ui ) =B(b, a) ui (1 + ui ) {1 + log(1 + 1/ui )}b+aπSB (ui ) =where a, b > 0 are hyperparameters and B(a, b) is the beta function. We discuss the rolesof the hyperparameters, a and b, of the proposed priors, and then we propose particular
choices of the hyperparameters. Under the proposed priors, we can construct an eﬃcient
sampling algorithm for posterior inference via Markov chain Monte Carlo (MCMC) method.
Furthermore, we show that the proposed priors have two desirable theoretical properties;(i) robust shrinkage rules for large observations, and(ii) Kullback-Leibler super-eﬃciency under sparsity.
The performance of the proposed method is demonstrated through numerical studies, andwe apply the method to two real datasets related to the average length of hospital stay for
COVID-19 in South Korea and variance estimation of gene expression data.

References[1] Carvalho, C. M., Polson, N. G. and Scott, J. G. (2010). The horseshoe estimator for sparse
signals. Biometrika, 97, 465–480.
[2] Datta, J. and Dunson, D. (2016). Bayesian inference on quasi-sparse count data.
Biometrika, 103(4), 971–983.
[3] Donoho, D. and Jin, J. (2006). Asymptotic minimaxity of false discovery rate thresholding
for sparse exponential data. The Annals of Statistics, 34, 2980–3018.
[4] Lu, M. and Stephens, M. (2016). Variance adaptive shrinkage (vash): ﬂexible empirical
Bayes estimation of variances. Bioinformatics, 32, 3428–3434.
[5] Polson, N. G. and Scott, J. G. (2010). Shrink globally, act locally: Sparse bayesian regularization and prediction. Bayesian Statistics 9, 501–538.

Minimum information dependence modelingfor mixed domain data
Keisuke Yano (the Institute of Statistical Mathematics),
Tomonari Sei (the University of Tokyo)
Abstract
In this talk, we introduce a joint statistical model for mixed-domain data that is
proposed by [1]. The proposed model contains multivariate Gaussian and log-linear
models. We show the existence and uniqueness of the proposed model under fairly
weak conditions. To estimate the dependence parameter in our model, we present
a conditional inference together with a sampling procedure and show it provides aconsistent estimator of the dependence parameter.

1

Minimum information dependence model

Let (Xi = 1, . . . , d be measure spaces and denote their product space by
Qi , F(Xi ), dxi ) for Q
the removal
X = di=1 Xi and dx = di=1 dxi . For index i, useQthe notation −i to indicate
Qof the i-th coordinate, e.g., x−i = (xj )j̸=i , X−i = j̸=i Xj , and dx−i = j̸=i dxj .
Let r1 (x1 ; ν), . . . , rd (xd ; ν) be statistical models of marginal densities on X1 , . . . , Xd ,respectively, where ν denotes parameters characterizing the marginal densities. We can assign, if necessary, independent parameters to each ri as ri (xi ; νi ) by setting ν = (ν1 , . . . , νd ).
We consider a class of probability density functions
! d
d
X
Y
p(x; θ, ν) = exp θ⊤ h(x) −ai (xi ; θ, ν) − ψ(θ, ν)ri (xi ; ν),
(1)
i=1

i=1

where θ ∈ RK is a K-dimensional parameter representing the dependence, and h : X → RK
is a given function. The functions ai (xi ; θ, ν) and ψ(θ, ν) are simultaneously determined byconstraintsZp(x; θ, ν)dx−i = ri (xi ; ν), i = 1, . . . , d, and
(2)
Z X
d

ai (xi ; θ, ν)p(x; θ, ν)dx = 0.

(3)

i=1

Note that the density (1) is reduced to the independent model

Qd

i=1 ri (xi ; ν) if θ = 0.

Definition 1. A statistical model (1) together with the constraints (2) and (3) is called aminimum information dependence model. The parameter θ is called the canonical parameter, ν is the marginal parameter, h(x) comprises the canonical statistics, ai (xi ; θ, ν)s are
the normalizing functions and ψ(θ, ν) is the potential function.
1

Figure 1 displays a two-dimensional histogram of samples from the minimum information dependence model for mixed variables (discrete and [0, 1]) with negative correlation,
which shows that the minimum information dependence model easily expresses the dependence between mixed-domain variables.
(a)

=0

y

0.0
9
8
7
6
5
4
3
2
1
0
0.0

0.2

0.2

0.4

0.4

0.6

x

0.6

(b)

= 100

0.0

0.2

0.4

0.0

0.2

0.4

0.8
9
8
7
6
5
4
3
2
1
0

0.8

x

0.6

0.8

0.6

0.8

Figure 1: Two-dimensional histograms of 10000 samples from the minimum information
dependence model with the Beta Beta(10, 10) and Poisson Po(3) marginals. The canonicalstatistic h(x, y) is given by h(x, y) = x/(y + 1). The joint histogram and marginal histograms are plotted. (a) Joint histogram with θ = 0. (b) Joint histogram with θ = 100.
Q
Let p0 (x) := di=1 ri (xi ; ν). We say that a function H ∈ L1 (p0 (x)dx) is feasible if there
exist measurable functions{ai (xi ) : i = 1, . . . , d} and a real number ψ ∈ R such that theP
H(x)− di=1 ai (xi )−ψfunction p(x) = ep0 (x) satisfiesZp(x)dx−i = ri (xi ; ν)for each i = 1, . . . , d andZ X
dai (xi )p(x)dx = 0.

i=1Theorem 1 (Theorem 1 of [1]). A function H ∈ L1 (p0 (x)dx) is feasible if there exist{bi ∈ L1 (ri (xi )dxi ) : i = 1, . . . , d} such thatZ
Pd
eH(x)− i=1 bi (xi ) p0 (x)dx < ∞.
(4)
Furthermore, if H is feasible, thenPd

i=1 ai (xi ) and ψ are unique.

References[1] T. Sei and K. Yano.
arXiv:2206.06792.

Minimum information dependence modeling,

2

2022.

高次元データ学習における特徴学習の優位性
Taiji Suzuki1,2
Joint work with Sho Okumoto , Jimmy Ba3 , Murat A. Erdogdu3 , Zhichao Wang4 ,Denny Wu3 , Greg Yang5
1

1

Graduate School of Information Science and Technology, the University of Tokyo,
2
RIKEN Center for Advanced Intelligence Project3
University of Toronto and Vector Institute, 4 University of California, San Diego,
5
Microsoft Research AI

In my presentation, I have presented two topics: (1) Learnability of convolutional neural networks for inﬁnite dimensional input, and (2) High-dimensional asymptotics of feature learning.
Each topic is based on the papers Okumoto & Suzuki (2022) and Ba et al. (2022) respectively.

1

Learnability of convolutional neural networks for inﬁnite
dimensional input

In the ﬁrst part, I have presented a recent result about the learning ability of convolutional neuralnetworks for inﬁnite dimensional input Okumoto & Suzuki (2022). Let λ be the uniform probability measure on ([0, 1], B([0, 1])) where B([0, 1]) is the Borel σ-ﬁeld on [0, 1], and let λ∞ be theproduct measure of λ on ([0, 1]∞ , B([0, 1]∞ )) where B([0, 1]∞ ) is the product σ-algebra generated
by the cylindric sets ∩j≤d {x ∈ [0, 1]∞ : xj ∈ Bj } for d = 1, 2, . . . and Bj ∈ B([0, 1]). Let PX be aprobability measure deﬁned on the measurable space ([0, 1]∞ , B([0, 1]∞ )) that is absolutely contindPX
uous to λ∞ and its Radon-Nikodym derivative satisﬁes ∥ dλ∞ ∥L∞ ([0,1]∞ ) < ∞. Then, suppose thato∞there exists a true function f : [0, 1] → R, and consider the following nonparametric regression
problem with an inﬁnite dimensional input:
Y = f o (X) + ξ,
where X is a random variable taking its value on [0, 1]∞ and obeys the distribution PX introduced
above, and ξ is a observation noise generated from N (0, σ 2 ) (a normal distribution with mean 0
and variance σ 2 > 0). Let P be the joint distribution of X and Y obeying the regression model.
Then, we analyzed the estimation accuracy of deep convolutional network to estimate f o from n
input-output observations (Xi , yi )ni=1 and investigated how the smoothness of the target function
f o aﬀects the convergence rate.
Among a wide range of success of deep learning, convolutional neural networks have been extensively utilized in several tasks such as speech recognition, image processing, and natural language
processing, which require inputs with large dimensions. Several studies have investigated function
estimation capability of deep learning, but most of them have assumed that the dimensionality of
the input is much smaller than the sample size. However, for typical data in applications such as
those handled by the convolutional neural networks described above, the dimensionality of inputs
is relatively high or even inﬁnite. In our analysis, we investigated the approximation and estimation errors of the (dilated) convolutional neural networks when the input is inﬁnite dimensional.
Although the approximation and estimation errors of neural networks are aﬀected by the curse ofdimensionality in the existing analyses for typical function spaces such as the Hölder and Besov
spaces, we show that, by considering anisotropic smoothness, they can alleviate exponential dependency on the dimensionality but they only depend on the smoothness of the target functions. Our
theoretical analysis supports the great practical success of convolutional networks. Furthermore,
we show that the dilated convolution is advantageous when the smoothness of the target function
has a sparse structure.
1

2High-dimensional asymptotics of feature learningt

In the second part, I have presented the asymptotic analysis of predictive accuracy of two layer
neural networks with feature learning developed in Ba et al. (2022).
We consider the training of a fully-connected two-layer neural network (NN) with N neurons,
1 ∑
1fNN (x) = √ai σ(⟨x, wi ⟩) = √ a⊤ σ(W ⊤ x),
N i=1
N
N

where x ∈ Rd , W ∈ Rd×N , a ∈ RN , σ is the nonlinear activation function applied entry-wise,
and the training objective is to minimize the empirical risk. Our analysis will be made in the
proportional asymptotic limit, i.e., the number of training data n, the input dimensionality d, andthe number of neurons N jointly tend to inﬁnity. Intuitively, this regime reﬂects the setting where
the network width and data size are comparable, which is consistent with practical choices of model
scaling.
We showed that the ﬁrst gradient update contains a rank-1 “spike,” which results in an alignment between the ﬁrst-layer weights and the linear component of the teacher model f ∗ . To characterize the impact of this alignment, we compute the prediction risk of ridge regression on the
conjugate kernel after one gradient step on W with learning rate η, when f ∗ is a single-index
model. We consider two scalings of the ﬁrst step learning rate η. For small η, we establish a Gaussian equivalence property for the trained feature map, and prove that the learned kernel improves
upon the initial random features model, but cannot defeat the best linear model on the input.
Whereas for suﬃciently large η, we prove that for certain f ∗ , the same ridge estimator on trainedfeatures can go beyond this “linear regime” and outperform a wide range of random features androtationally invariant kernels. Our results demonstrate that even one gradient step can lead to aconsiderable advantage over random features, and highlight the role of learning rate scaling in theinitial phase of training.

References
J. Ba, M. A. Erdogdu, T. Suzuki, Z. Wang, D. Wu, and G. Yang. High-dimensional asymptotics
of feature learning: How one gradient step improves the representation. In Advances in Neural
Information Processing Systems, volume 35, 2022. to appear.
S. Okumoto and T. Suzuki. Learnability of convolutional neural networks for inﬁnite dimensional
input via mixed and anisotropic smoothness. In International Conference on Learning Representations, 2022.

2

量子乱数プロトコルを用いた量子計算機の安定性評価
鹿野 豊
群馬大学大学院理工学府
Institute for Quantum Studies, Chapman University
email: yshikano@gunma-u.ac.jp
概 要

All computing devices, including quantum computers, must exhibit that for a given input, an output is producedin accordance with the program. The outputs generated by quantum computers that fulﬁll these requirements are nottemporally correlated, however. In a quantum computing device comprising solid-state qubits such as superconducting
qubits, any operation to rest the qubits to their initial state faces a practical problem. On the implementation of
the scalable quantum computers, the health check (or stability check) algorithms are needed. We propose that thequantum random number generation is one of the candidates of the health check algorithms in any quantum computing
devices.

1

はじめに量子計算機に関して、量子アルゴリズムと呼ばれる計算アルゴリズムの根本的な変更により、これまで知られていた従来型の計算量とは別の量子計算量という概念を創出した。特筆すべきは、現代暗号の計算量的安全性の基盤として用いられている「素因数分解のアルゴリズム」は NP (Non-deterministic Polynomial time) のクラスであることは知られているが、量子コンピュータによって誤り確率が高々1/3 の多項式時間で解ける決定問題の複雑性クラスである BQP (Bounded-error

Quantum Polynomial time) で解けるという Simon-Shor のアルゴリズムというものが知られている。そのため、量子計算機を物理系として実現するための開発のレシピを必死に模索してきており、現在は小規模ではあるが、クラウド上でユーザーが自由にアルゴリズムを実行できる量子計算機も実現している。一方で、量子状態をそのまま情報のソースとして用いるためにノイズに対して従来型のデジタル計算機より耐性がないことが知られている。このため、一つ一つの計算素子の単位である量子ビットの正確な評価を行い、物理的な原因を追跡した上で量子ビットの改善をこれまで行ってきた。しかし、このようなボトムアップでのやり方は量子ビット数が増えれば増えるほど、評価だけに時間がかかってしまうため、ユーザーが計算機として使える時間が減ってしまうということが知られている。そのため、簡易に量子計算機の評価・診断ができるプロトコルの開発が必要であるということが必要であることは必然であろう。
そこで、本論稿では、量子計算機を計算機システムの一つだと捉えると、インプットおよびアウトプットにはビット列を用いたブラックボックスモデルとシステム的な観点で考えることが出来る。この際、デジタル計算機と量子計算機の本質的な差は、インプットに対してアウトプットが確定的に出てくるデジタル計算機に対して 1 、インプットに対してアウトプッ
トが確率的に出てくるため、何度も同じアルゴリズムを実行し、その結果から統計的な判断を加えなければならないのが量子計算機の特徴である。すると、量子計算機の特徴を調べるためには、本質的に確率的な出力を行うシステムであるということを統計的に判定する問題に落とし込むことが出来る。そのため、本論稿では、量子計算機が確率的出力を本質的にするシステムであることをチェックするためにシードのない乱数生成器 [1]2 であるということを検証する方法を検討する。

2

量子乱数プロトコル
我々の提案する量子乱数プロトコルは非常に単純である。量子状態を用意する空間として量子力学の公理より 2 次元の複素ヒルベルト空間 H := C2 を考える。そして、初期の量子状態として( )1|0⟩ :=∈ C2
0(1)

1 当然、デジタル計算機にも必須なノイズがあるため、誤り訂正符号という技術を実装させ、確定的にアウトプットのビット列を出力させなければならない。
2 シードのある乱数生成器のことの代表例はデジタル計算機の中で実装される疑似乱数生成器である。物理乱数生成器の多くは初期条件の違いに対して出力結果が鋭敏であることを利用しているため、真性乱数生成器とは呼ばれていない。

1

を用意する。そして、ここにアダマールゲート H ∈ U (2) という状態を作用させると
( )
1
1 1
|ψ⟩ := H|0⟩ = √ (|0⟩ + |1⟩) = √
2
2 1

(2)

となる。量子力学の公理より測定結果に関しては、量子状態に対して自己共役作用素 A をスペクトル分解する。今回の場
合、A = σz ∈ U (2) というパウリ行列を採用すると、(
σz = +1|0⟩⟨0| + (−1)|1⟩⟨1| = +1
†

1 0
0 0

)

(
+ (−1)

0 0
0 1)=: +1P̂0 + (−1)P̂1(3)†となる。ここで、⟨0| := (|0⟩) , ⟨1| := (|1⟩) という Dirac 記法を採用している。これは、射影作用素 P̂0 , P̂1 に対して、それぞれに対応する固有値が測定結果として得られる。その固有値が得られる確率はボルンの確率規則より

Pr (∗) := ∥P̂∗ |ψ⟩∥2 , (∗ ∈ {0, 1})(4)となることから、今回の例においては1(5)2となる。そして、測定結果を得た後に初期状態にリセットさせる。こうすることで何度でもアルゴリズムが実行できるよう
Pr (0) = Pr (1) =になる。これは 1 量子ビットに対して、全て同じ操作を行っていることからシードのない乱数生成器と捉えることが出来る。
本論稿では、これを用意されたクラウド量子計算機の全ての量子ビットに対して、同じ操作を「同時に」3 。次に、このアルゴリズムを 5 日間実行させ、その結果として量子ビットが安定的な確率的な出力を得ているかということを単純に実行回数に対して「1」の出た積算回数をグラフ化することにより直線になっていなければ不安定に動作したということが分かった [3]。そして、初期状態へのリセットが不十分な場合、出力されたビット列に対して時間相関が検出されてしまう可能性があるため、チェックを行ったところ、時間相関がランダムに行っている可能性が示唆された [4]。このような検定を行う
ことで量子計算機の安定動作に関する指標を構築できる可能性があるということを示した。

3

おわりに
量子計算機をシステムとして捉えることで、量子乱数プロトコルが量子計算機の簡易診断プロトコルになるということを指摘した。本論稿で提案したプロトコルをクラウド量子計算機で実行したところ、全て安定的に動作されていないということが統計的な処理により明らかになった。
今後、NIST Test Suites を拡張させ、バイアスのある確率出力がある場合での乱数検定手法を開発させ、更には、安定性動作指標として本論稿では実行回数に対して「1」の出た積算回数をグラフ化することにより直線になっていなければ不安定に動作したということを判定基準としたが、その他にもランダム化試験のような統計的手法を適用することが重要である。統計学的基準を利用した量子計算機の指標づくりは、これまでも行われてきたが、今後、更に重要であると予想するため理論的な側面のみならず、実際のデータ処理をした際のソフトウェアの整備まで含めて開発していかねばならない宿命で
ある。

参考文献
[1] Y. Shikano, AIP Conference Proceedings 2286, 040004 (2020).
[2] K. Tamura and Y. Shikano, TUCS Lecture Notes 30, 13 - 25 (2019); Cryptology ePrint Archive, Paper 2020/078.
[3] K. Tamura and Y. Shikano, International Symposium on Mathematics, Quantum Theory, and Cryptography, Mathematics for Industry, vol 33 (Springer, Singapore, 2021) 17– 37.
[4] Y. Shikano, K. Tamura, and R. Raymond, EPTCS 315, 18– 25 (2020).

3 ハードウェアの制約上、同時に本当に実行されているわけではなさそうであるが、残念ながら量子ビット数を増やして動作させるという方向性に昨
今は焦点が当てられており、どれくらい正確に量子操作が実行できるかという点においてはあまり改善が見られていない。特に長期的に安定的に動作できるシステムではないということは強調しておきたい。更には、2022 年 10 月現在であっても全ての乱数検定をパスするクラウド量子計算機は存在していない。

2

Bootstrap for Selecting a subset which contains all populations
better than a standard
Jun-ichiro Fukuchi, Gakushuin University, Tokyo, Japan
1 Framework
When samples are obtained from k populations with the diﬀerent means, statistician is often interested
in ﬁnding populations whose means are larger than that of the standard population. For the case of normal
populations, Gupta and Sobel (1958) formulated subset selection approach. In this approach, the goal is toselect a subset which contains all populations better than a standard. At this presentation, we only describe
the method for unknown standard case. The method for known standard is similar and it is omitted. Papers
related to this methodology are Paulson (1952), Dunnet (1955), Tong (1969), Huang, Panchapakesan and Tseng(1984) among others.
A related problem is selecting a subset which contains the population with the largest mean. The basic
framework is built by Gupta (1965). Swanepoel (1983) investigated the bootstrap method for subset selectionof the largest mean. Fukuchi (2020) proved consistency and the second-order correctness of the bootstrap under
the assumption that population distributions belong to the same location family. This type of homogeneity
assumption is needed for the selection of largest mean problem since the statistic used for selection rule depends
on the distribution function associated with the largest mean asymmetrically. Subset selection problem of better
than standard is simpler in some sense and as a result the bootstrap is shown in this paper to be successful
without homogeneity assumption on the distributions even when sample sizes are unequal.

2 Selection rule and the inﬁmum of the probability of the correct selection
Let Π1 , . . . , Πk be k populations with means µ1 , . . . , µk , and the common variances σ 2 . Let Π0 be thestandard population. Denote the distribution function of Πi by Fi , i = 0, 1, . . . , k. It is assumed that thefunctional forms of Fi are unknown. Let µ[1] ≤ µ[2] ≤ · · · ≤ µ[k] be ordered means. Let {Xij : j = 1, . . . , ni }is a sample from the population Πi (i = 0, 1, . . . , k). Let X i be the sample mean of {Xij , j = 1, 2, . . . , ni },i = 0, 1, . . . , k and S 2 be the pooled variance.

2.1 Equal sample sizes case
Consider the following selection rule proposed by Gupta and Sobel (1957).
Selection rule R1 : Retain the population Πi in the selected subset if and only if
c1
Xi ≥ X0 − √ Sn(1)where c1 is the smallest constant such that inf Ω P (CS|R1 ) = p for prespeciﬁed p ∈{(µ0 , µ1 , . . . , µk ) : µi ∈ R, i = 0, . . . , k}. We make following assumption.

( −1 )
k , 1 . Let Ω :=Assumption 1. For each i = 0, . . . , k, Fi belongs to a location family, that is Fi (x) = Gi (x − µi ) for some
distribution function Gi .
Theorem 1. Under Assumption 1,(
inf P (CS|R1 ) =PΩ
Ω

(
max

1≤i≤k

X 0 − µ0X −µ√ − i√ i
S/ n
S/ n

)

)
≤ c1

2.2 Unequal sample sizes caseIn this section, we consider the unequal sample size case. Consider the following selection rule.
Selection rule R2 : Retain the population Πi in the selected subset if and only if
√
1
1
+ ,
X i ≥ X 0 − c2 Sn0niwhere c2 is the smallest constant such that inf Ω P (CS|R2 ) = p for prespeciﬁed p ∈ (k −1 , 1). See Bechhofer,
Santner and Goldsman (1995), p130 for this selection rule for normal populations.
Theorem 2. Under Assumption 1,

inf P (CS|R2 ) =PΩ  max
Ω

1≤i≤k

(X 0 − µ0 ) − (X i − µi )√≤ c2 
S n10 + n1i
1

(2)

)−1
(
)( √
Let T2,n := max1≤i≤k (X 0 − µ0 ) − (X i − µi ) S n10 + n1i. Then c2 is the p-th quantile of T2,n where p
is prespeciﬁed number.

3 Estimation of c2 by bootstrap
Estimating c2 = c2 (p) amounts to estimating the distribution function of T2,n . Since Fi ’s are unknown, theconstant c = cn (p) needs to be estimated from the sample. Let n = (n0 , n1 , . . . , nk ) and∗
Xn = {Xij : j = 1, . . . , ni , i = 0, 1, . . . , k}. Let Xn∗ = {Xij: j = 1, . . . , ni , i = 0, 1, . . . , k, } be the bootstrap∗∗
sample. Let P denote the conditional distribution of Xij, j = 1, . . . , ni , i = 1, . . . , k given Xn . Based on the
)( √
)−1
( ∗
∗
∗
is
bootstrap sample Xn∗ , the bootstrap version T2,n:= max1≤i≤k (X 0 − X 0 ) − (X i − X i ) S ∗ n10 + n1i
computed.
The( ∗) bootstrap estimator ĉ2 (p) of c2 (p) is deﬁned by the p-th quantile of the bootstrap distribution
P ∗ T2,n≤ x . Theorem 3 states that bootstrap method described above is consistent. Unlike the bootstrap
for subset selection of the largest mean, the bootstrap for subset selection of better than standard is consistentunder rather weak assumptions on the population distributions and the rates how sample sizes increase.
Theorem 3. Let Xi1 , . . . , Xini be i.i.d. random variables with E(Xij ) = µi , V(Xij ) = σ 2 , i = 0, . . . , k.
Assume that0k(i) {X0j }nj=1, . . . , {Xkj }nj=1
are independent.
(ii) There exist di ∈ [0, ∞], i = 1, . . . , k such that ni /n0 → di .
Then, as min ni → ∞0≤i≤k∗sup |P ∗ (T2,n≤ x) − P (T2,n ≤ x)| → 0, a.s.

(3)

x∈R

and |ĉ2 (p) − c2 (p)| → 0, a.s. for any p ∈ (1/k, 1).
When sample sizes are equal and population distributions are absolutely continuous, by Theorem 3 of Fukuchi(2020), bootstrap estimation of T2,n is second-order correct meaning that bootstrap approximation is better
than normal approximation.

4 Concluding remarks
In this presentation, properties of bootstrap when it is used to estimate the constant in subset selection rule of
selecting better than standard are investigated. We assumed the common unknown variance. Consistency andsecond-order correctness for the proposed bootstrap method are proved. Results of simulation study is given
at the presentation.

References[1] Bechhofer, R. E., Santner, T. J. and Goldsman, D. M. (1995). Design and Analysis of Experiments for
Statistical Selection, Screening, and Multiple Comparisons. Wiley.
[2] Fukuchi, J. (2020). A note on bootstrap for Gupta’s subset selection procedure. Sankhya A, volume 82, p.
96–114 (2020).
[3] Gupta, S. S. (1965). On Some Multiple Decision (Selection and Ranking) Rules. Technometrics, 7, No. 2,
225-245[4] Gupta, S. S. and Hsiao, P. (1983). Empirical bayes rules for selecting good populations. Journal of
Statistical Planning and Inference. Volume 8, Issue 1, September 1983, Pages 87-101[5] Gupta, S. S. and Sobel, M. (1958). On selecting a subset which contains all populations better than astandard. Annals of Mathematical Statistics Vol. 29, No. 1 (Mar., 1958), pp. 235-244.
[6] Huang, D.Y., Panchapakesan, S. and Tseng, S. T. (1984). Some locally optimal subset selection rules forcomparison with a control. Journal of Statistical Planning and Inference. Volume 9, Issue 1, January 1984,
Pages 63-72[7] Paulson, E. (1952) On the Comparison of Several Experimental Categories with a Control. The Annals of
Mathematical Statistics, Vol. 23, No. 2 (Jun., 1952), pp. 239-246.
[8] Tong, Y. L. (1969). On partitioning a set of normal populations by their locations with respect to acontrol. Ann. Math. Statist. 40, 1300-1324.
[9] Swanepoel, J. W. H. (1983). Bootstrap selection procedures based on robust estimators. Comm. Statist.
Theory Methods, 12(18), 2059-2083.

2