Պ‫ݚ‬අγϯϙδ΢Ϝʮ౷‫ܭ‬తϞσϦϯάͱ‫ࢉܭ‬ΞϧΰϦζϜͷ਺ཧͱల։ʯ
‫ج‬൫‫ڀݚ‬ʢAʣ 15H01678ʮେ‫ن‬໛ෳࡶσʔλͷཧ࿦ͱํ๏࿦ͷ૯߹త‫ڀݚ‬ʯ‫ڀݚ‬୅දऀɿ ੨ౢ ੣ (ஜ೾େֶ)։࠵੹೚ऀɿ ۚ৿ ‫ܟ‬จ (໊‫ݹ‬԰େֶ)೔࣌ɿ 2017 ೥ 2 ݄ 18 ೔ (౔)ʙ19 ೔ (೔)৔ॴɿ ໊‫ݹ‬԰େֶ ৘ใՊֶ‫ڀݚ‬Պ౩ϓϩάϥϜ2/18 (౔)

09:55 - 10:00 Φʔϓχϯά੨ౢ੣ (ஜ೾େ)10:00 - 11:40 (25 ෼ ×4)
஛೭಺ߴࢤ (͸ͩͯ͜ະདྷେ/ཧ‫ ݚ‬AIP) ม‫ܗ‬ϒϨάϚϯٖ‫ͦͱ཭ڑ‬ͷपล໼ా࿨ળ (ஜ೾େ) ߴ࣍‫ݻݩ‬༗஋ɾ‫ݻ‬༗ϕΫτϧͷҰகਪఆྔʹ͍ͭͯ‫ޖ‬Ғ (໊େ) ૊߹ͤ໰୊ʹର͢Δ min-max regret ‫ج‬४ͷϩόετ࠷దԽࠓઘҸ૱ (౦େ) ࡞༻ૉਪఆɿؔ਺σʔλճ‫ؼ‬ͷ౷‫ܭ‬తੑ࣭

13:10-14:50 (25 ෼ ×4)
٢ᖒਅଠ࿠ (‫఼ޚ‬৔‫ )ݚૅج‬Ϧοδճ‫ؼ‬ͷ૒ରද‫ݱ‬খ࡚හ‫( ׮‬εςϥϦϯΫ (‫ ))ג‬ճ‫ؼ‬෼ੳͱਲ਼‫ܭ‬ըङ໺मฏ (౦େ) Aﬃne ‫ʹࣜܗ‬ΑΔϗϩϊϛοΫޯ഑๏ͷ‫ࠩޡ‬ධՁ໦ଜ‫۝( ࣇܓ‬େ) ม਺બ୒ʹର͢Δࠞ߹੔਺ඇઢ‫ܭܗ‬ը๏

15:10-16:20(25 ෼ ×1 & 45 ෼ (ಛผߨԋ)×1)

ࠇՏఱ (౦େ) ଟ֯‫ॆܗ‬రߏ଄Λ࣋ͭը૾σʔλ΁ͷ࠷దάϥϑຒΊࠐΈख๏ͷ։ൃྛਖ਼ਓ (໊େ/NUS) ಛผߨԋɿInformation Geometry Approach to Parameter Estimation in Hidden Markov Models

16:40-18:20 (25 ෼ ×4)
தࢁ༏‫( ޗ‬ஜ೾େ) ߴ࣍‫ݩ‬খඪຊʹ͓͚ΔόΠΞεิਖ਼ SVM
ߐ‫ޱ‬ਅಁ (౷਺‫ )ݚ‬Generalized Boltzmann machines and activation functions‫۽‬୩࿱ (ਆಸ઒େ) ύϥϝʔλసҠֶशͷཧ࿦ղੳদౡ৻ (౦େ) େ‫ن‬໛ͳઢ‫ܗ‬༧ଌ‫ث‬ͷͨΊͷ ඇಉ‫ظ‬ಛ௃நग़εΩʔϜ

2/19 (೔)
10:00 - 11:40 (25 ෼ ×4)
ฏඌক߶ (Ѫ஌‫ཱݝ‬େ) ߦྻࣜ఺աఔͷ४ϞϯςΧϧϩੵ෼΁ͷԠ༻ᔉۚ๕ (ઍ༿େ) Cell Regression and Reference Priors
ᖒਖ਼‫( ݑ‬ਆ‫ށ‬େ) Euclidean Design Theory੺ؒཅೋ (౦๺େ) TBA

13:10-14:50 (25 ෼ ×4)
দҪ޹ଠ (໊େ) ೖΕࢠ‫߹ࠞܕ‬Ϟσϧʹ‫ ͮ͘ج‬cancer outlier profile ͷਪఆͱ͕Μ਍அ΁ͷԠ༻ʹ͍ͭͯ
ҏ౻௚‫( ل‬౦େ) ࣮༻తͳՃ଎ۙ઀ޯ഑๏ͷ࣮૷ͱ 2 ஋൑ผϞσϧ΁ͷԠ༻ຊ୩ल‫޻໊( ݎ‬େ) ۭؒͷ௿ϥϯΫੑͱฏ‫׈‬ੑΛߟྀͨ͠ ϑʔϦΤ܎਺࠷దԽʹΑΔ MR
௒ղ૾ҏ౻৳Ұ (౦େ) େࣗ༝౓‫ܥ‬ͷσʔλಉԽͷͨΊͷ 2nd-order adjoint ๏Λ༻͍ͨߴ଎ෆ֬
࣮ੑධՁ๏

15:10-16:50 (25 ෼ ×4)
ਗ਼ஐ໵ (౦େ) ࠲ඪ͝ͱͷม਺ม‫ʹ׵‬ΑͬͯಘΒΕΔ Stein ‫ܕ‬ͷ෼෍୞໦޹ଠ࿠ (த෦େ) An operational characterization of the notion of probability by algorithmic randomness and its applicationsத઒݈࣏ (௕ԬٕՊେ) ௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍ͷࣹӨΞϧΰϦζϜʹΑΔ୳ࡧʹ͍ͭͯ౉ลҰൕ (๛‫ٕڮ‬Պେ) Ϩʔτ࿪Έཧ࿦ͱҰൠԽࣄ‫ޙ‬෼෍

16:50 - 17:00 Ϋϩʔδϯά変形ブレグマン擬距離とその周辺
竹之内高志
はこだて未来大学, 理化学研究所 革新知能統合研究センター
概要
離散空間上の確率モデルは様々な現象・データを柔軟に表現可能な有用なツールであり, データから適切にパラメーターを決定することで精度の高い推論が可能となる. その一方で, 離散空間上の確率モデルは正規化項の計算に指数オーダーの計算量を必要とすることがあるため, パラメーターの推定量を構成することが難しいことが多い. 本研究では, 変形ブレグマンダイバージェンスと e-混合モデルを用いることで正規化項の計算をすることなく構成可能な推定量を提案し, その性質について議論する.

1

導入

X を d 次元の離散空間 X 上の確率変数ベクトルとし, ⟨f ⟩ =る. 本稿では確率モデルq̄θ (x) =!

qθ (x), qθ (x) = exp(ψθ (x)), Zθ = ⟨qθ ⟩Zθ

x∈X f (x) とす(1)に着目し, パラメーター θ を推定することを目的とする. ただし, qθ (x) は非正規化モデル, Zθ = ⟨qθ ⟩ はモデル q̄θ が確率測度であることを要請する正規化項である. データセット D = {xi }ni=1 に対応する経験分布を p̃(x) とする.
"
nx
x ∈ Z,
p̃(x) = n0otherwise,
ただし Z = {x|x ∈ D}, nx はデータセット D に含まれる x の個数とする. 最尤推定量は経験分布 p̃(x) と確率モデル q̄θ (x) 間の KL ダイバージェンス最小化として定式化することができるが, 正規化項の計算に由来する計算量の問題
から推定量の構成が困難であることが多い. 本稿では変形ブレグマンダイバー
ジェンスと e-混合モデルを用いて, 正規化項の計算を行わずに構築可能な推定量を提案する.

2

提案法

2 つの正値測度 f, g に対して, 以下の変形ブレグマンダイバージェンスを考える. ただし U は凸関数, f は単調関数とする.
#$
D(p, q; U, f ) = U (f (q)) − U (f (p)) − U ′ (f (p)) {f (q) − f (p)} .
(2)
1

D(p, q; U, f ) ≥ 0 が常に成立し, D(p, q; U, f ) = 0 が成立するのは p = q のときのみである (Murata et al., 2004). 経験分布 p̃(x) と確率モデル q̄θ (x) の e混合モデルを以下のように定義する (Amari and Nagaoka, 2000) .
" % &αnxqθ (x)1−α x ∈ Z,p̃α (x)q̄θ (x)1−αp̃α (x)qθ (x)1−α
n
#
$
#
$
r̃α,θ (x) ==∝.
p̃α q̄θ1−αp̃α qθ1−α0otherwise,
ただし α(̸= 0, 1) を定数とし, r̃0,θ (x) = q̄θ (x), r̃1,θ (x) = p̃(x) が成立する. e混合モデルは正規化項を計算することなく構成可能であることに注意する.
本稿では, これらを用いて, 以下のような推定量を考える.

θ̂ U,f = argmin D(r̃α,θ , r̃α′ ,θ ; U, f ).

(3)

θe-混合モデルが正規化項の計算をすることなく構成可能であるため, (3) も正
規化項の計算を行うことなく構成可能な推定量となる. この推定量に対して,
以下が成立する.
命題 1 推定量 θ̂ U,f はフィッシャー一致性を持つ. すなわち, データを生成す
る真の分布が p(x) = q̄θ0 (x) であり, rα,θ =q̄θα (x)q̄θ (x)1−α0!
"とするとき, 以下q̄θα q̄θ1−α0が成立する.

θ 0 = argmin D(rα,θ , rα′ ,θ ; U, f ).

(4)

θ

補題 1 推定量 θ̂ U,f の漸近分布は漸近的に以下の多変量正規分布に従う.
√−1−1n(θ̂ U,f − θ 0 ) ∼ N (0, HU,f,θ
J
HU,f,θ
).
(5)
0 U,f,θ 0
0
'
(
#
$
T
ただし HU,f,θ0 = ξU,f (q̄θ0 )q̄θ0 (ψθ′ 0 − µ0 )(ψθ′ 0 − µ0 )T , JU,f,θ0 = q̄θ0 ζU,f,θ0 ζU,f,θ
,
0
#
$
′
′′
′
2
µ0 = q̄θ0 ψθ0 , ξU,f (z) = U (f (z))f (z) z,
ζU,f,θ (x) = ξU,f (q̄θ (x))(ψθ′ (x) − µ0 ) − ⟨q̄θ ξU,f (q̄θ )(ψθ′ − µ0 )⟩ である.
定理 1 関数 f がU ′′ (f (z))f ′ (z)2 z = 1

(6)

−1
−1
を満たす時, 推定量 θ̂ U,f の漸近分散 HU,f,θ
J
HU,f,θ
はフィッシャー情報
0 U,f,θ 00量行列の逆行列となり, 漸近有効となる.

References
S. Amari and H. Nagaoka. Methods of Information Geometry, volume 191
of Translations of Mathematical Monographs. Oxford University Press,
2000.
G. E. Hinton and T. J. Sejnowski. Learning and relearning in Boltzmann
machines. MIT Press, Cambridge, Mass, 1:282–317, 1986.
N. Murata, T. Takenouchi, T. Kanamori, and S. Eguchi. Information geometry of U -boost and bregman divergence. Neural Computation, 16(7):1437–1481, 2004.
T. J. Sejnowski. Higher-order Boltzmann machines. In American Instituteof Physics Conference Series, 151(7):398–403, 1986.
2

高次元固有値・固有ベクトルの一致推定量について
矢田
青嶋

和善
誠

(筑波大数理物質)(筑波大数理物質)

1. はじめに情報化の進展に伴い，高次元データの統計解析が，ますます重要になってきている．
2000 年以降，確率論と理論物理の方面から，ランダム行列の理論に基づく幾つかの重要な結果がもたらされた．Johnstone (2001, AS) 等は，標本固有値の漸近分布を導出した．しかしながら，そこでは，データの次元数 d と標本数 n が n/d → c > 0 を満たす場
合を考え，高次元において標本数は次元数と同程度を仮定した．例えば，次元数は優に 10,000 を超えるが標本数は高々100 程度といった高次元小標本においては，標本数を次元数と同程度には仮定できない．それゆえ，n が d に依存しないような設定で，もしくは，n = n(d) であっても n/d → 0 となる設定で，高次元漸近理論を展開する必要がある．Yata and Aoshima [2] は，高次元小標本における PCA の性質を研究し，PCAが一致性をもつための標本数 n の d に関するオーダー条件を導き，高次元小標本において PCA が不適解を起こすことを示した．この問題を解決する策として，Yata andAoshima [3] は，高次元小標本データ空間の幾何学的表現を研究し，それに基づいて “ノ
イズ掃き出し法”とよばれる方法論を考案した．一方で，Yata and Aoshima [4] は，高次元大標本も含む一般的な高次元データに対して，power spiked モデルと呼ばれる固有値モデルを考案し，高次元データに対する新しい PCA を構築した．最近，Aoshima
and Yata [1] は，ノイズ掃き出し法による固有ベクトルの推定量を用いることで新たな高次元二標本検定法を考案した．
本講演では，高次元固有ベクトルの一致性について論じ，閾値を用いてノイズ掃き出し法による固有ベクトルの推定量を補正することで，緩い仮定のもと高次元固有ベ
クトルの一致性を与える新たな方法論を提案した．

2. 高次元固有ベクトルの一致性

平均に d 次のベクトル µ，共分散行列に d 次の半正定値行列 Σ をもつ母集団を考える．母集団から n (≥ 3) 個の d 次データベクトル x1 , ..., xn を無作為に抽出する．
Σ の固有値を λ1 ≥ · · · ≥ λd (≥ 0) とし，適当な直交行列 H = [h1 , ..., hd ] で Σ を
Σ = HΛH T , Λ = diag(λ1 , ..., λd ) と分解する．標本共分散行列 S のスペクトル分解を!
T
S = di=1 λ̂i ĥi ĥi とする．最近，Yata and Aoshima [4] は，power spiked モデルとよばれる固有値モデルを考案し，高次元データに対する新しい PCA を研究した. いま，!
!pTTΣ(1) = m
j=1 λj hj hj ，Σ(2) =j=m+1 λj hj hj とおき，Σ = Σ(1) + Σ(2) という分解を考える．そのとき，次の条件を満たすような λ1 ≥ · · · ≥ λd を power spiked モデルと定義する．
m
λm に対して，limd→∞ tr(Σk(2))/λkmm = 0 なる (有界な) ある自然数 km が存在する. (1)いま，δj = λ−1j tr(Σ(2) )/(n − 1)，j = 1, ..., m とおく．power spiked モデル (1) のもと，次の定理を得る．
定理 1 ([4]).

各 j = 1, ..., m について，適当な正則条件と条件

(C-i) tr(Σ2(2) )2 /(nλ4j ) = o(1)のもと，d, n → ∞ のとき次が成り立つ．

λ̂j= 1 + δj + op (1) and hTj ĥj = (1 + δj )−1/2 + op (1).
λj定理 1 より，適当な正則条件と (C-i) のもと次を得る．

||ĥj − hj ||2 = 2{1 − (1 + δj )−1/2 } + op (1).

(2)

ここで，|| · || はユークリッドノルムを表す．一方で，Yata and Aoshima [3] は，高次元
小標本データ空間の幾何学的表現を研究し，それに基づいて “ノイズ掃き出し法”とよばれる方法論を考案し，次のような固有値の推定量を提案した．
!
tr(S) − ji=1 λ̂iλ̃j = λ̂j −(j = 1, ..., n − 2).
(3)
n−1−j

さらに，Σ の固有ベクトルについて，ノイズ掃き出し法による推定を考える．推定量 (3)に基づいて，Σ の固有ベクトル hj を h̃j = (λ̂j /λ̃j )1/2 ĥj で推定する．そのとき，power
spiked モデル (1) のもと，次の定理を得る．
定理 2 ([4]). 各 j = 1, ..., m について，適当な正則条件と (C-i) のもと，d, n → ∞ のとき次が成り立つ．
λ̃j= 1 + op (1) and hTj h̃j = 1 + op (1).
λj
それゆえ，hj の内積に関する一致性をもつ．
ここで，||h̃j ||2 = λ̂j /λ̃j ≥ 1 であることに注意し，定理 1 と 2 より，適当な正則条件と (C-i) のもと次を得る．
||h̃j − hj ||2 = δj + op (1).
(4)
すなわち，(2) と (4) より，lim inf d,n→∞ δj > 0 のとき，ĥj と h̃j はノルムに関する一致性をもたない．
本講演では，閾値を用いて h̃j を補正することで，δj → ∞ のもとでも，||h́j − hj ||2 = op (1)が成り立つような新たな固有ベクトルの推定量 h́j を提案し，理論的かつ数値的に既存の推定量と比較した．

参考文献[1] Aoshima, M., Yata, K. (2016). Two-sample tests for high-dimension, strongly spiked
eigenvalue models. Statist. Sinica, in press (arXiv:1602.02491).
[2] Yata, K., Aoshima, M. (2009). PCA consistency for non-Gaussian data in high dimension,
low sample size context. Commun. Statist. Theory Methods 38 2634-2652.
[3] Yata, K., Aoshima, M. (2012). Eﬀective PCA for high-dimension, low-sample-size datawith noise reduction via geometric representations. J. Multivariate Anal. 105 193-215.
[4] Yata, K., Aoshima, M. (2013). PCA consistency for the power spiked model in highdimensional settings. J. Multivariate Anal. 122 334-354.

૊߹ͤ໰୊ʹର͢Δ min-max regret ‫ج‬४ͷϩόετ࠷దԽ
Nagoya University Wei WU
Nagoya University Mutsunori YAGIURA

Several optimization problems arising in real world applications do not have accurate estimates of the problem parameters when the optimization decision is taken. Stochastic programming and robust optimization are two common approaches for the solution of optimization
problems under uncertainty. The min-max and min-max regret criteria are two of the typicalapproaches for robust optimization. The min-max criterion aims at obtaining a solution withthe best worst-case value across all scenarios. The regret is defined as the diﬀerence between
the actual cost and the optimal cost that would have been obtained if a diﬀerent solution had
been chosen. The min-max regret approach is to minimize the worst-case regret. This criterion
is not as pessimistic as the min-max approach.
For the max-min regret criterion, we consider the knapsack problem (KP) under discrete
profits. For the discrete scenario case, we assume that the scenario set is described explicitly. The
discrete max-min knapsack problem (MM-KP) is known to be strongly NP-hard for unboundedscenario set [3]. However, the MM-KP is solvable by a pseudo-polynomial time algorithm when
the size of scenario set is bounded by a constant. We examine this pseudo-polynomial time
method based on dynamic programming. We also examine a branch-and-cut algorithm based amixed integer programming (MIP) model. We propose a heuristic algorithm for the MM-KP thatsolves the underlying KP to optimality under a fixed scenario. For the average-profit scenario,
we show that the optimal value of this fixed-scenario KP under the average-profit scenario isa valid upper bound. We further propose an iterative method to improve the performance ofthe fixed scenario heuristic. In each iteration, we generate a new scenario based on solutions
obtained by then.
For the min-max criterion, we consider the generalized assignment problem (GAP) and themultidimensional knapsack problem (MKP) under interval costs. The classical GAP is a strongly
NP-hard combinatorial optimization problem [6] having many applications (see [1, 4], and [5]).
The classical MKP is a strongly NP-hard combinatorial optimization problem [2] and has been
widely studied over many decades due to both theoretical interests and its broad applications inseveral engineering fields, such as cargo loading, cutting stock, bin-packing, financial and othermanagement issues [7].
The interval min-max regret generalized assignment problem (MMR-GAP) is a generalizationof the GAP to the case in which the cost coeﬃcients are uncertain. In real life applications, thecosts are often aﬀected by many factors, and they can be unknown at the optimization stage. We
assume that every cost coeﬃcient can take any value in a corresponding given interval, regardlessof the values taken by the other cost coeﬃcients. The problem requires to find a robust solutionthat minimizes the maximum regret. We prove that the decision version of MMR-GAP isΣp2 -complete. We propose a heuristic algorithm for the MMR-GAP that solves the underlying
GAP to optimality under a fixed scenario. We consider three scenarios (lowest cost, highest
cost, and median cost), and we show that the median cost scenario leads to a solution of the
MMR-GAP whose objective function value is within twice the optimal value. We also propose1

a dual substitution heuristic based on a MIP model obtained by replacing some constraints
with the dual of their continuous relaxation. We also propose exact algorithmic approaches
that iteratively solve the problem by only including a subset of scenarios. The first approachis based on logic-based Benders decomposition: it solves a MIP with incomplete scenarios anditeratively supplements the scenarios corresponding to violated constraints. We then introduce
a basic branch-and-cut algorithm and enhance it through: (i) Lagrangian relaxations, to provide
tighter lower bounds than those produced by the linear programming relaxation; (ii) an eﬃcientvariable fixing technique; (iii) a two-direction dynamic programming approach to eﬃciently solve
the Lagrangian subproblems. We compare the introduced algorithms through computational
experiments on diﬀerent benchmarks.
For the interval min-max regret multidimensional knapsack problem (MMR-MKP), we propose a new heuristic framework, which we call the iterated dual substitution (IDS) algorithm.
The IDS iteratively generates linear constraints (rows) based on a mixed integer programming
model. Computational experiments on a wide set of benchmark instances are carried out, andthe proposed iterated dual substitution algorithm performs best on all of the tested instances.

References[1] M.L. Fisher, R. Jaikumar, “A generalized assignment heuristic for vehicle routing,” Networks, 11 (1981) 109–124.
[2] M. Garey, D. Johnson, Computers and Intractability: A Guide to the Theory of NPcompleteness, W.H. Freeman, San Francisco, 1979.
[3] P. Kouvelis and G. Yu, Robust Discrete Optimization and Its Applications, Kluwer, Academic Publishers, Dordrecht, 1997.
[4] S. Martello, P. Toth, Knapsack Problems: Algorithms and Computer Implementations, John
Wiley & Sons, Chichester, New York, 1990.
[5] K.S. Ruland, “A model for aeromedical routing and scheduling,” International Transactions
in Operational Research, 6 (1999) 57–73.
[6] S. Sahni, T. Gonzalez, “P-complete approximation problems,” Journal of the Association for
Computing Machinery, 23 (1976) 555–565.
[7] M. Varnamkhasti, “Overview of the algorithms for solving the multidimensional knapsackproblem,” Advanced Studies in Biology, 4 (2012) 37–47.

2

Operator Estimation : Analysis for Functional Regression
with Functional Input and Output
Masaaki Imaizumi (UT)(joint work with Kengo Kato (UT))

In this presentation, we investigate a regression problem where both covariate
and response variables are random functions. Let the covariate X and the response
Y be L2 (I)-valued random variables with I = [0, 1], and T : L2 (I) ! L2 (I) be aregression operator. Then, we consider that the random functions are generated
from the following regression model
Y (t) = T (X)(t) + ✏(t), t 2 I,

(1)

where ✏ is a Gaussian noise process in L2 (I).
Data representation by functions is used for analyzing data which are observed ona large number of grids in its domain. The fields of statistical methods for analyzing
such data is called functional data analysis, and it is summarized in Ramsay andSilverman (2005).
The functional regression with functional covariates and functional responses isinvestigated by many studies, especially the linear regression case. Cuevas et al.
(2002), Yao et al. (2005), Crambes and Mas (2013), Hörmann and Kidziński (2015)
and more studies propose their estimators and clarify their theoretical
R properties.
Since it is possible to represent the linear regression operator by a form b(·, t)x(t)dt
using a bivariate function b(s, t), some of the methods for the linear regression are
conducted by estimating b(s, t).
The nonlinear regression with functional covariate and functional response is a
developing problem, since there is no standard way to represent the nonlinear operator T . Such the problem is considered by Bosq and Delecroix (1985), and various
methods are proposed. The Nadaraya-Watson method with the kernel function isinvestigated by Ferraty et al. (2011), Lian (2011), Ferraty et al. (2012), and others.
Another method for the problem is a functional reproducing kernel Hilbert space(fRKHS) method which is studied by Preda (2007), Lian (2007), and Kadri et al.
(2015).
In this presentation, we propose an estimator for the linear and nonlinear functional regression problem, and clarify some regularity conditions and properties ofthe estimators. Then, we derive the convergence rate of the estimator which is characterized by smoothness of the variables and the operator, where the smoothness
represents a speed of delay of coefficients measured by the basis function decomposition. About the convergence rate of the estimator for the linear regression, we

1

find that the rate is independent from the smoothness of the operator in the direction of responses, and the rate is identical to the convergence rate with the linear
regression with a scalar output case provided by Hall and Horowitz (2007). For the
nonlinear case, we also find the similar properties of the convergence rate with someassumptions.

References
Bosq, D. and Delecroix, M. (1985) Nonparametric prediction of a hilbert spacevalued random variable, Stochastic processes and their applications, 19, 271–280.
Crambes, C. and Mas, A. (2013) Asymptotics of prediction in functional linear
regression with functional outputs, Bernoulli, 19, 2627–2651.
Cuevas, A., Febrero, M. and Fraiman, R. (2002) Linear functional regression: thecase of fixed design and functional response, Canadian Journal of Statistics, 30,
285–300.
Ferraty, F., Laksaci, A., Tadj, A. and Vieu, P. (2011) Kernel regression with functional response, Electronic Journal of Statistics, 5, 159–171.
Ferraty, F., Van Keilegom, I. and Vieu, P. (2012) Regression when both response
and predictor are functions, Journal of Multivariate Analysis, 109, 10–28.
Hall, P. and Horowitz, J. L. (2007) Methodology and convergence rates for functional
linear regression, The Annals of Statistics, 35, 70–91.
Hörmann, S. and Kidziński, L. (2015) A note on estimation in hilbertian linear
models, Scandinavian journal of statistics, 42, 43–62.
Kadri, H., Duflos, E., Preux, P., Canu, S., Rakotomamonjy, A. and Audi↵ren, J.
(2015) Operator-valued kernels for learning from functional response data, Journalof Machine Learning Research, 16, 1–54.
Lian, H. (2007) Nonlinear functional models for functional responses in reproducing
kernel hilbert spaces, Canadian Journal of Statistics, 35, 597–606.
Lian, H. (2011) Convergence of functional k-nearest neighbor regression estimate
with functional responses, Electronic Journal of Statistics, 5, 31–40.
Preda, C. (2007) Regression models for functional data by reproducing kernel hilbert
spaces methods, Journal of Statistical Planning and Inference, 137, 829–840.
Ramsay, J. and Silverman, B. (2005) Functional Data Analysis, 2ne Edition,
Springer.
Yao, F., Müller, H.-G. and Wang, J.-L. (2005) Functional linear regression analysis
for longitudinal data, The Annals of Statistics, 33, 2873–2903.

2Ϧοδճ‫ؼ‬ͷ૒ରද‫ݱ‬٢ᖒ ਅଠ࿠

(‫఼ޚ‬৔‫ૅج‬Պֶ‫ڀݚ‬ձ)∗֓ཁϦοδճ‫ؼ‬ͷ̍ͭʹ Hodrick-Prescott ϑΟϧλ (1997) ͱ‫ݺ‬͹ΕΔ΋ͷ͕͋Δ
ʢҎԼɺHP ϑΟϧλͱ‫Ϳݺ‬ʣɻ͜ͷϑΟϧλ͸ Leser(1961) ͕‫͞ͱݯى‬Ε͍ͯ
ΔɻHP ϑΟϧλ͸ɺ‫؍‬ଌ࣌‫ྻܥ‬ϕΫτϧͱτϨϯυ੒෼ϕΫτϧͷࠩͷ̎৐ϊϧϜʹɺࠩ෼ߦྻ͕τϨϯυ੒෼ϕΫτϧʹ࡞༻ͨ͠ϕΫτϧͷ̎৐ϊϧϜΛεΧϥഒͯ͠Ճ͑ͨ΋ͷͱͳ͍ͬͯΔɻYamada(2015) ͸ɺHP ϑΟϧλͷτϨϯυ੒෼ϕΫτϧʹɺࠩ෼ߦྻ͔Βߏ੒͞ΕΔม‫ͨ͠༻࡞͕׵‬଒ʹରͯ͠ɺ۩ମతͳ 2 छྨͷ࡞༻ʹ͍ͭͯߟ࡯͍ͯ͠Δɻ͜ͷ 2 छྨͷ࡞༻͸ɺτϨϯυ੒෼ϕΫτϧʹ࡞༻͢Δࠩ෼ߦྻͷۭ֩ؒͱ૾ۭ͔ؒΒ۩ମతʹߏ੒͞Εͨ΋ͷͰ͋ΔɻຊདྷɺτϨϯυ੒෼ϕΫτϧ΁ͷม‫ʹͯ΂͢༻࡞׵‬ର͠ɺϦοδճ‫ؼ‬ͷද‫ݱ‬Λ‫ݕ‬౼͢΂͖Ͱ͋Δͱߟ͑Δɻͦ͜ͰɺຊใࠂͰ͸ࠩ෼ߦྻͷۭ֩ؒͱ૾ۭؒͱΛཅʹදࣔ͢Δ͜ͱͰɺHP ϑΟϧλͷҰൠද‫ݱ‬‫ͦͼٴ‬ͷ૒ରද‫ݱ‬Λߟ࡯͢Δɻ͜ΕΒ̎ͭͷҰൠද‫ݱ‬ͷϖΞ͔Β௨ৗͷ૒ର‫ز‬Կߏ଄ʢ৘ใ‫ز‬ԿʣΛߟ࡯ͨ͠ɻͨͩ͠ɺ௨ৗͷϦοδճ‫ؼ‬ͷਖ਼ଇԽύϥϝʔλ͸ඇෛΛલఏͱ͍ͯ͠Δ͕ɺҰൠද‫ͦͱݱ‬ͷ૒ରද‫͍͓ͯʹݱ‬͸ɺਖ਼ଇԽύϥϝʔλ͕ඇෛͱ͸‫ݶ‬Βͳ͍ͱ͜Ζ͕௨ৗͷਖ਼ଇԽͱ૬ҧ͍ͯ͠Δɻ
ਖ਼ଇԽύϥϝʔλ͕ෛͷ஋΋औΓ͏Δ͜ͱ͸ɺ‫ز‬Կֶతʹ͸ෆࣗવͳ͜ͱͰ
͸ͳ͍ɻਖ਼ଇԽύϥϝʔλͱࠩ෼ߦྻͷಛҟ஋ͷେখؔ܎͕ʢҰൠද‫ͼٴݱ‬ͦͷ૒ରද‫ݱ‬ͷϔογΞϯ͕ਖ਼ఆɺෆఆɺෛఆͷ̏ͭͷ৔߹ʣɺHP ϑΟϧλͷੑೳʹ༩͑ΔӨ‫ڹ‬Λ਺ཧతʹߟ࡯͢Δ͜ͱ͸ɺࠓ‫ޙ‬ͷ՝୊ͱ͢Δɻࢀߟจ‫ݙ‬[1] C.E.V.Leser, A simple method of trend construction, Journal of the Royal StatisticalSociety, Series B (Methodological), 23, (1961), 91-107.
[2] R.J.Hodrick and E.C. Prescott, Postwar U.S. business cycles: An empirical investigation,
Journal of Money, Credit and Banking, 29(1),(1997), 1-16
[3] H.Yamada, Ridge regression representations of the generalized Hodrick-Prescott filter, J.
Japan Statist.Soc., Vol.45, No.2 (2015), 121-128.

2016 ೥౓ Պ‫ݚ‬අγϯϙδ΢Ϝɿ౷‫ܭ‬తϞσϦϯάͱ‫ࢉܭ‬ΞϧΰϦζϜͷ਺ཧͱల։ʢԙɿ໊‫ݹ‬԰େֶ৘ใֶ‫ڀݚ‬Պ, 2016 ೥ 2 ݄ 18 ೔ʙ19 ೔ʣ∗e-mail: yzw2003@mail.goo.ne.jp

報告書 (回帰分析と錐計画)
小崎敏寛 (Toshihiro Kosaki)∗
ステラリンク株式会社 (Stera Link, Co., Ltd.)

1

課題
回帰分析は，一つの被説明変数の変動をいくつかの説明変数によって説明する手法である．その中で，あてはめからの誤差からなるベクトルの 2 ノルムの 2 乗を最小にするようにパラメータを決める手法が最小二乗法である．最小化する基準として，1 ノルムや ∞ ノルムを考えることも古くから行われてきた[2, 8]．

2

解決手法
近年，錐計画問題に注目が集まっている．錐計画問題の一つのクラスとして，二次錐計画問題 [1, 6] がある．最近，この問題の 2 ノルムを一般のノルムにしたノルム錐計画問題 [4, 5] を提案した．本研究では，ノルムをさらに一般化したゲージ関数 [3, 7] を考える．ゲージ関数からなる錐を持つ問題をゲージ
錐計画問題と呼ぶ．最小化する基準として，2-ノルムでなく，ゲージ関数を考える．すると問題はゲー
ジ錐計画問題になる．

3

結果
ゲージ錐計画問題の双対問題を考えることで，弱双対定理がなりたつことを示した．

4

まとめと今後の課題
回帰分析の計算法にゲージ錐計画が適用できることを示した．この問題に対して，弱双対定理がなりたつことを示した．
今後の課題として次のようなものがある．ゲージ錐計画を解くアルゴリズムを考え，実際に数値実験を行う．そして得られた結果の解釈を行う．また，双対ギャップが存在しない条件を調べる．

参考文献[1] F. Alizadeh and D. Goldfarb, Second-Order Cone Programming, Mathematical Programming, 95,
3-51, 2003.
∗ toshihirokosaki@gmail.com

1[2] T. S. Arthanari and Y. Dodge, Mathematical Programming in Statistics, John Wiley & Sons, Inc.,
1981.
[3] H-H Chao and L. Vandenberghe, Semideﬁnite Representations of Gauge Functions for Structured
Low-Rank Matrix Decomposition, arXiv, 2016.
[4] 小崎敏寛, ノルム錐計画問題の双対性, 京都大学数理解析研究所講究録 1931：最適化アルゴリズムの進展：理論・応用・実装, 89-93, 2015.
[5] 小崎敏寛, 一般のノルム錐計画問題の弱双対定理, 統計数理研究所共同研究リポート 369, 最適化：モデリングとアルゴリズム 28, 75-78, 2016．

[6] M. Lobo, L. Vandenberghe, S. Boyd, and H. Lebret, Applications of Second-Order Cone Programming, Linear Algebra and its Applications, 284, 193-228, 1998.
[7] D. G. Luenberger, Optimization by Vector Space Methods, John Wiley & Sons, Inc., 1969.
[8] G. A. Watson, Approximation in Normed Linear Spaces, Journal of Computational and Applied
Mathematics, 121, 1-36, 2000.

2

Aﬃne 形式によるホロノミック勾配法の誤差評価
東京大・情報理工 狩野 修平
東京大・情報理工 清智也

Nakayama et al. (2011) によって提案されたホロノミック勾配法は，計算対象の関数をそれが満たす微分方程式系に変換して数値計算を行うという手法である．ホロノミック
勾配法を適用する場合，微分方程式を解くために適当な初期値を与える必要がある．既存研究の多くではパラメータ付き積分があるパラメータにおいて陽に書ける場合はその値を用い，それが不可能な場合は級数展開を用いて初期値を設定している (Sei and Kume(2015) 等)．しかし，級数展開による初期値設定は真の初期値を計算できているわけではなく，実際に計算した積分値がどれぐらい信頼できるのかということはあまり考察されてこなかった．その問題に対処するために，本発表では初期値が区間で与えられたもとでのPfaﬃan 系の数値計算について考える．
数値計算したいホロノミック関数を f (θ) = f (θ1 , . . . , θp ), θ ∈ Θ ⊂ Rp とし，f (θ) を含む q 次のベクトルを m(θ) と書く．i = 1, . . . , p について，m(θ) に∂i m(θ) = Pi (θ)m(θ),

!

∂
∂i =
∂θi

"

(1)

を満たすような Pfaﬃan 系 Pi : Θ → Rq×q が存在しているものとする．この時，興味のあるパラメータを θ(1)，それと異なるパラメータを θ(0) と書き，それらをつなぐ curveを θ(t), t ∈ [0, 1] と表す．この時式 (1) は∂t m(θ(t)) = Pt∗ m(θ(t)) = g(t, m),Pt∗ :=p#(∂t θi (t))Pi (θ(t)),

∂t =

i=1

∂
∂t

(2)

のような常微分方程式系に書き換えられる．式 (2) で m(θ(0)) が与えられた時に m(θ(1))
を計算することを考える．ここで，与えられた初期値 m̂(θ(0)) については，|m̂(θ(0)) − m(θ(0))| ≤ C, C = (c1 , . . . , cq )⊤ ∈ Rq

(3)のような評価が成立していると仮定する．式 (2)(3) で与えられた微分方程式に対するアプローチとしては，Taylor model による古典的な方法 (Lohner (1987) 等) と，既存の Runge-Kutta 法のスキームを利用して解の包含を生成する手法が提案されている(Dit Sandretto and Chapoutot (2016) 等)．本研究では，ホロノミック勾配法の既存研究において幅広く Runge-Kutta 法が用いられていることを考慮し，後者の枠組みを用いて微分方程式を数値的に解き，既存研究と比較することを考える．

1

式 (2)(3) で表される微分方程式の数値解を計算するには，区間演算と呼ばれる変数の上界と下界を与えて区間を表す演算法が用いることができる．しかし，区間演算は変数間の相関性を反映できないという欠点があるため，本研究では Aﬃne 形式 (Stolfi and De

Figueiredo (2003)) と呼ばれる変数の存在範囲を多項式の形で表す手法を使って計算を行う．区間演算と Aﬃne 形式は，演算の前後に対して包含性が保存されるという性質があるため，初期値を区間の形で与え，その中に真の値が含まれている限り，計算結果として出力される区間に真の解が含まれているとみなすことができる．具体的には，t = tj における Aﬃne 形式の m̃j とあるステップ幅 h について，次のように tj+1 = tj + h における

Aﬃne 形式 m̃j+1 を計算する．
˜m̃j+1 = ỹn+1 + LTE([tj , tj+1 ], m([tj , tj+1 ]))
ただし，ỹn+1 は

ỹn+1 = m̃n + hs#bi k i ,ki = g(tn + ci h, m̃n + h

i=1

s
#

aij kj )j=1のように，通常の Runge-Kutta 法で計算を行い，LTE([tj , tj+1 ], m) は各ステップにおける打ち切り誤差を示す．発表では実際に用いる計算アルゴリズムの詳細と数値実験の結果について説明する．

参考文献
Dit Sandretto, J. A. and Chapoutot, A. (2016). Validated Explicit and Implicit RungeKutta Methods, Reliable Computing, Vol. 22, pp. 78–103.
Lohner, R. J. (1987). Enclosing the Solutions of Ordinary Initial and Boundary Value
Problems, in Computer Arithmetic: Scientific Computation and Programming Languages, Stuttgart: Wiley-Teubner, pp. 255–286.
Nakayama, H., Nishiyama, K., Noro, M., Ohara, K., Sei, T., Takayama, N., andTakemura, A. (2011). Holonomic gradient descent and its application to the FisherBingham integral, Advances in Applied Mathematics, Vol. 47, No. 3, pp. 639–658.
Sei, T. and Kume, A. (2015). Calculating the normalising constant of the Bingham
distribution on the sphere using the holonomic gradient method, Statistics andComputing, Vol. 25, No. 2, pp. 321–332.
Stolfi, J. and De Figueiredo, L. (2003). An Introduction to Aﬃne Arithmetic, TEMA
- Tendências em Matemática Aplicada e Computacional, Vol. 4, No. 3, pp. 297–312.
2

変数選択に対する混合整数非線形計画法
九州大学大学院数理学府数理学専攻 木村 圭児
九州大学マス・フォア・インダストリ研究所 脇 隼人

1.

変数選択と AIC 最小化

統計学におけるモデル推定では, 情報量規準を用いてモデルに採用するパラメータを選択することがあり, これは変数選択と呼ばれる. 赤池情報量
規準 (Akaike’s information criterion: AIC) などの情報量規準が用いられ, AIC が最小であるようなパラメータの組合せを求めることで, データとの当てはまりの良さを損なわないように予測精度の良いモデルを推定することができる. AIC ができるだけ小さい値をとる変数選択の一般的な手法の一つとして, ステップワイズ法が知られている.
ステップワイズ法は, R などの既存の統計ソフト
ウェアに実装されていて, 計算が高速で広く用いられている. しかし, ステップワイズ法は局所探索をしているとみなせるので, 推定されたモデルのAIC が最小であるとは限らない.

2.

既存手法と提案する手法

ロジスティック回帰モデルに対して, 線形近似を用いた混合整数線形計画法による変数選択が提案されている [3]. この手法によって得られる解は,AIC が最小であると保証されないが, ステップワ
イズ法に比べ質の良い解である. 一方, 線形回帰モ
デルに対する変数選択の手法として, 混合整数二次錐計画法を用いた手法 [2] や混合整数二次計画法を用いた手法 [1] が提案されている.
本研究は, 最適化を基にした変数選択に対して,
混合整数非線形計画問題として定式化し, 定式化された問題を効率良く解くための様々な工夫を提案している. 提案する手法を実装するために, 分
枝限定法のフレームワークを提供している SCIP
(Solving Constraint Integer Program [4]) を用いる. SCIP は自由度が高いソフトウェアであり, 解法を細かく制御するプログラムを実装することができる. そのため, 本研究は定式化した問題をソルバーのみで解くというわけでは無く, 効率良く工夫を提案し, 実装まで実現している. 線形回帰における AIC 最小化に対して, 提案する手法がどの程度の規模の問題まで求解できるか紹介する.

3.

MINLP として定式化

最適化を用いる変数選択では, 一般に, 目的関数
は “与えられたデータとモデルの誤差” と “説明変数の数” の二つの項から形成される. モデルに現れるパラメータを β = (β1 , . . . , βp )T ∈ Rp とす
る. 変数選択では, 説明変数の候補が選択されなかった場合, 対応するパラメータ βj は 0 である.
変数選択のための最適化問題は, 次の混合整数非線形計画問題として定式化できる.
!
min f (β) + λ
zj
β,z

j∈Ip

s.t. zj = 0 ⇒ βj = 0 (j ∈ Ip )βj ∈ R, zj ∈ {0, 1} (j ∈ Ip )(1)

ただし, λ は正の定数で, Ip = {1, . . . , p} である.
f (β) が与えられたデータとモデルの誤差を表す関数ならば, 問題 (1) の目的関数はモデルを評価する"p
ことができる. 目的関数の第二項 λ j=1 zj は説
明変数の数のペナルティとして機能する.

4.

下界値と上界値の計算

分枝限定法を用いた問題 (1) を効率良く手法を提案する. 分枝限定法は, 問題 (1) の部分問題の最小値の下界値と問題 (1) の最小値の上界値が必要である. 分枝限定法のフレームワークを提供している SCIP[4] に下界値と上界値の計算を実装する
ことで, 問題 (1) を解くことができる.
分枝限定法における分枝操作によって, 部分問題が生成される際, zj (j ∈ Ip ) は 1 あるいは 0 に固定される. 生成された部分問題に対して, zj (j ∈ Ip )の添字に関して次の集合を定義する.

Z1 = {j ∈ Ip : zj は 1 に固定されている },Z0 = {j ∈ Ip : zj は 0 に固定されている },Z = {j ∈ Ip : zj はまだ固定されていない }.

部分集合 Z1 , Z0 , Z ⊆ Ip に対して, 問題 (1) の部

分問題 Q(Z1 , Z0 , Z) は次のように表される.
!
min f (β) + λ
zj
β,z

5.

j∈Ip

s.t. zj = 0 ⇒ βj = 0, zj ∈ {0, 1} (j ∈ Z)
βj = 0, zj = 0 (j ∈ Z0 )zj = 1 (j ∈ Z1 ), βj ∈ R (j ∈ Ip )(2)

部分問題 Q(Z1 , Z0 , Z) の zj ∈ {0, 1} (j ∈ Z) の整数性を緩和した緩和問題は次のように表される.
!
min f (β) + λ
zj
β,z

j∈Ip

s.t. zj = 0 ⇒ βj = 0, 0 ≤ zj ≤ 1 (j ∈ Z)
βj = 0, zj = 0 (j ∈ Z0 )zj = 1 (j ∈ Z1 ), βj ∈ R (j ∈ Ip )(3)

問題 (3) の最小値は, 部分問題 Q(Z1 , Z0 , Z) の最小値の下界値である. 次に問題 (3) から制約 zj = 0 ⇒βj = 0, 0 ≤ zj ≤ 1 (j ∈ Z) と zj (j ∈ {1, . . . , p})を取り除いた以下の問題について考える.

min f (β) + λ#(Z1 )βs.t. βj = 0 (j ∈ Z0 ), βj ∈ R (j ∈ Ip )(4)

問題 (3) の任意の実行可能解 (β̄, z̄) に対して,""j∈Ip z̄j =j∈Z z̄j + λ#(Z1 ) より, 次の不等式が成立する.
!
f (β̄) + λz̄j ≥ f (β̄) + λ#(Z1 ).

効率良く解くための工夫

本研究では, 問題 (1) を効率良く解くために以下の工夫を提案している.
(I) 親問題の緩和問題を用いた下界値の更新(II) 解の傾向を利用した計算コストの削減(III) 分枝変数の決め方(IV) ステップワイズ法を基にした上界値の更新
分枝限定法のフレームワークを提供している
SCIP[4] にこれらを実装している.

6.

数値実験

線形回帰における AIC 最小化に対して, 既存の手法 (MISOCP[2] と MIQP[1]) と提案する手法(MINLP) の実験結果を紹介する. 数値実験には,[5] で公開されているデータを使用する. 計算時間
は制限時間を 5000 秒とする. 5000 秒で解けない
場合は, >5000 と記し, AIC は得られた上界値を
記す.
name
p 手法
AIC time(sec)

sfC

26

fires

63

crime

100

j∈Ip

β̄ は (4) の実行可能解である. よって, 問題 (4) の最小値は部分問題 Q(Z1 , Z0 , Z) の最小値の下界値である. したがって, 本研究では, 問題 (4) を部分問題 Q(Z1 , Z0 , Z) の下界値を求めるための緩和問題として扱う.
線形回帰における AIC 最小化の場合, 問題 (4)は制約無し凸二次計画問題となる. よって, 線形方程式を解くことで問題 (4) の最小解を求めることができる. ロジスティック回帰における AIC 最小化の場合, 問題 (4) は制約無し凸計画問題となる.
よって, ニュートン法などの既存のアルゴリズムを適用することで問題 (4) の最小解を求めることができる.
分枝限定法では, 問題 (1) の最小値の上界値も必要である. 問題 (4) の最小解から問題 (1) の実行可能解を生成することができる. 生成した解を用いて, 問題 (1) の最小値の上界値を得る.

MINLP
MISOCP[2]
MIQP[1]MINLP
MISOCP[2]
MIQP[1]MINLP
MISOCP[2]
MIQP[1]2816.32816.32816.31429.61431.61435.13410.33690.73646.4

10.5
489.8
26.49
>5000
>5000
>5000
>5000
>5000
>5000

参考文献
[1] D. Bertsimas, A. King and R. Mazumder, “Best
Subset Selection via a Modern Optimization
Lenz” , Ann. Stat., 44, 2, 813 – 852, 2016.
[2] R. Miyashiro and Y. Takano, “Mixed Integer
second-order cone programming formulations forvariable selection,” Eur. J. Oper. Res., 247, 721–
731, 2015.
[3] T. Sato, Y. Takano, R. Miyashiro and A. Yoshise,“Feature Subset Selection for Logistic Regression
via Mixed Integer Optimization,” Computational
Optimization and Applications, 2015.
[4] SCIP: Soving Constraint Integer Programs, http://scip.zib.de/[5] UCI Machine Learning Repository, http://archive.ics.uci.edu/ml/

多角形充填構造を持つ画像データへの最適グラフ埋め込み手法の開発
黒河天 ∗伊藤伸一 †

長尾大道 ∗†

糟谷正 ‡

井上純哉 ‡

1 はじめに細胞や組織の画像から，それらの境界線を抽出する問題は，多くの分野で重要である．例えば，材料科学分
野では，材料内部における結晶粒成長の写真から粒の境界を抽出し，データ同化を用いたパラメータ推定や状態予測に生かすことができる．結晶粒の成長モデルに適用可能なデータ同化手法として [1] がある．領域分
割 [2] やクラスタリング [3] といった，各ピクセルにラベルを与えるアプローチは，境界線そのものを求める
ことができない．一方，エッジ検出や直線検出は，物理的に解釈可能な線群を与えるとは限らない．
頂 点 集 合 V ，辺 集 合 E の 定 め る 無 向 グ ラ フ G =(V, E ) と ２ 次 元 標 本 点 上 の ス カ ラ ー 場

f = (fm,n )m=1,...,M,n=1,...,N （ 例 え ば ，写 真 の 画 素 値 な ど ）に 対 し て ，f を最 も よ く 説 明 す る G の 平
面埋め込みを求める問題をグラフ当てはめ (graph fitting)，略して graphit と呼ぶことにする．我々はgraphit の観点から，多角形充填構造を持つ画像の境界線を求める手法を提案する．

2 提案手法
２次元標本点上のスカラー場 f を連続領域上に拡張したスカラー場 Ψ を定義する．頂点 V = {v1 , . . . , vI }の座標 (x1 , y1 ), . . . , (xI , yI ) が定める G の平面埋め込みを考え，埋め込みに沿った線積分

J (x1 , y1 , . . . , xI , yI ) =!

{vi ,vj }∈E

"

Ψ dl

(1)

Γij

を導入する．ただし，Γij は (xi , yi ) から (xj , yj ) への線分を表す．このとき，最適化問題

min J (x1 , y1 , . . . , xI , yI ) s.t. (x1 , y1 , . . . , xI , yI ) ∈ S(2)の解を，求める graphit とする．ただし，S は適当な実行可能領域を表す．また，Ψ ではなく Ψ 2 の線積分を最小化する方法を，２乗 graphit という．
線分 Γij を関数

# $
# $
xi
xj
Φij (θ) = (1 − θ)+θ, θ ∈ [0, 1]yiyj∗ 東京大学大学院情報理工学系研究科† 東京大学地震研究所‡ 東京大学大学院工学系研究科

1(3)

でパラメータ付けるとき，graphit における J の勾配は陽に書けて

∂J
=
∂xi
∂J
=
∂yi

!

j:{vi ,vj }∈E!

j:{vi ,vj }∈E((

xi − xj
L2ij
yi − yj
L2ij"

Ψ dl + Lij

Γij"Ψ dl + Lij

Γij

" 1
0

" 1
0

◦ Φ )(θ) dθ),( ∂Ψ∂x(1 − θ)ij(1 − θ)ij( ∂Ψ∂y ◦ Φ )(θ) dθ)

(4)
(5)

となる．Lij は Γij の長さとする．したがって graphit は逐次２次計画法 [4] や L-BFGS-B [5] で計算できる．

Ψ の定義には幾つかの方法が考えられる．
1. 以下の Ψ : [0, M ) × [0, N ) → R の定める graphit を，素朴な graphit という：Ψ |[m−1,m)×[n−1,n) (x, y ) := fm,n for m = 1, . . . , M, n = 1, . . . , N.

(6)

2. 以下の Ψ : [0, M − 1) × [0, N − 1) → R の定める graphit を，１次スプラインに基づく graphit という：

Ψ |[m−1,m)×[n−1,n) (x, y ) :=(m − x)(n − y )fm,n + (m − x)(y − n + 1)fm,n+1+ (x − m + 1)(n − y )fm+1,n + (x − m + 1)(y − n + 1)fm+1,n+1 (7)for m = 1, . . . , M − 1, n = 1, . . . , N − 1.

それぞれの定める graphit 及び２乗 graphit を用いて，多角的充填構造を持つ人工画像の境界を抽出する数値実験を行った．この結果，１次スプラインに基づく２乗 graphit が，最も遠い初期頂点配置から真の頂点配置を探索する傾向にあった．

3 おわりに
本研究では graphit に基づく画像からの境界抽出の手法を提案した．提案手法は必ず多角形充填構造を持つ解を与えるが，解は頂点の初期配置に強く依存するため，手動で引いた境界線を補正するに止まる．今後の課題として，良い初期配置を自動的に構成する手法の開発が挙げられる．

参考文献[1] S. Ito, H. Nagao, A. Yamanaka, Y. Tsukada, T. Koyama, M. Kano, and J. Inoue. Data assimilation
for massive autonomous systems based on a second-order adjoint method. Physical Review E, Vol. 94,
No. 043307, 2016.
[2] N. R. Pal and S. K. Pal. A review on image segmentation techniques. Pattern Recognition, Vol. 26,
pp. 1277–1294, 1993.
[3] A. K. Jain, M. N. Murty, and P. J. Flynn. Data clustering: a review. ACM Computing Surveys,
Vol. 31, pp. 264–323, 1999.
[4] P. T. Boggs and J. W. Tolle. Sequential quadratic programming. Acta Numerica, Vol. 4, pp. 1–51,
1995.
[5] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained
optimization. SIAM Journal on Scientific Computing, Vol. 16, pp. 1190–1208, 1995.
2

Information Geometry Approach to
Parameter Estimation in Hidden
Markov Models
Masahito Hayashi1
Graduate School of Mathematics, Nagoya University, Japan,
and Centre for Quantum Technologies, National University of Singapore,
Singapore.
1

Abstract: We consider the estimation of hidden Markovian process byusing information geometry with respect to transition matrices. We consider
the case when we use only the histogram of k-memory data. Firstly, we
focus on a partial observation model with Markovian process and we showthat the asymptotic estimation error of this model is given as the inverseof projective Fisher information of transition matrices. Next, we apply thisresult to the estimation of hidden Markovian process. For this purpose, we
define an exponential family of Y-valued transition matrices. We carefully
equivalence problem for hidden Markovian process on the tangent space.
Then, we propose an novel method to estimate hidden Markovian process.

Information geometry established by Amari and Nagaoka [1] is a very powerful
method for statistical inference. Recently, the paper [4] applied this approach toestimation of Markovian process. In the paper [4], they employed informationgeometry of transition matrices given by Nakagawa and Kanaya [2] and Nagaoka[3]. Since this geometric structure depends only on the transition matrices, itdoes not change as the number n of observation increases while the geometry
based on the probability distribution changes according to the increase of thenumber n. In particular, the paper [4] introduced the curved exponential familyof transition matrices, and derived the Cramér-Rao inequality for the family,
which shows the optimality of the inverse of the transition matrix version ofFisher information matrix.
On the other hand, some of preceding studies [5, 6] of hidden Markov process
employed information geometry. However, they studied em-algorithm based on
the geometry of probability distributions, which changes according to the increase of the number n. So, the estimation process becomes complicated whenn is large. Hence, they could not evaluate the asymptotic behavior of the estimation error.
In this paper, we apply the information geometry of transition matrices toestimation of hidden Markovian process. Since we need to estimate the hidden
structure from the observed value, we apply the em algorithm based on thegeometry of transition matrices. That is, for this purpose, we formulate a partial
observation model of Markovian process and the em algorithm based on the1

M. Hayashi/Information Geometry Approach in Hidden Markov Models

2

geometry of transition matrices for this model. Then, using the transition matrix
version of the projective Fisher information, we evaluate the asymptotic error
in the Markovian case under a certain regularity condition.
However, we have another difficulty for the hidden Markovian process. Thereis ambiguity for the transition matrix to express the hidden Markovian process.
That is, there is a possibility that two di↵erent transition matrices express
the same hidden Markovian process. This problem is called the equivalence
problem and is solved by Ito, Kobayashi, and Amari [7]. However, to discuss theestimation of the hidden Markovian process, we need to consider this problem
in the tangent space because the asymptotic error is characterized by the local
geometrical structure.
In this paper, for this purpose, we establish the local equivalence relation forthe hidden Markovian process. When we apply the above em algorithm to thehidden Markovian process, for the regularity condition, we need to guaranteethat the tangent space is non-degenerate with respect to the local equivalence
condition.
Acknowledgment
The authors are very grateful to Professor Takafumi Kanamori and Professor
Vincent Y. F. Tan for helpful discussions and comments. The works reported
here were supported in part by the JSPS Grant-in-Aid for Scientific Research(B) No. 16KT0017, the Okawa Research Grant and Kayamori Foundation ofInformational Science Advancement.
References[1] S. Amari and H. Nagaoka, Methods of Information Geometry. Oxford University Press (2000).
[2] K. Nakagawa and F. Kanaya, “On the converse theorem in statistical hypothesis testing for markov chains,” IEEE Trans. Inform. Theory, Vol. 39,
No. 2, 629-633 (1993).
[3] H. Nagaoka, “The exponential family of Markov chains and its informationgeometry” Proceedings of The 28th Symposium on Information Theory and
Its Applications (SITA2005), Okinawa, Japan, Nov. 20-23, (2005).
[4] M. Hayashi and S. Watanabe, “Information Geometry Approach to Parameter Estimation in Markov Chains,” Annals of Statistics, Volume 44,
Number 4, 1495-1535 (2016).
[5] S. Amari, “Information geometry of the EM and em algorithms for neuralnetworks,” Neural Networks, Vol. 8, No. 9, 13791408 (1995).
[6] Y. Fujimoto and N. Murata, “A modified EM algorithm for mixture models based on Bregman divergence,” Annals of the Institute of Statistical
Mathematics, Vol. 59, No. 1, 325 (2007).
[7] H. Ito, S. -I. Amari, and K. Kobayash, “Identifiability of Hidden Markov Information Sources and Their Minimum Degrees of Freedom,” IEEE Trans.
Inform. Theory, Vol. 38, No. 2, 324-333, (1992).

高次元小標本におけるバイアス補正 SVM
筑波大学・数理物質科学

中山優吾

筑波大学・数理物質系

矢田和善

筑波大学・数理物質系

青嶋誠

1 はじめに
本講演では，高次元小標本データにおける判別分析を考えた．母集団が 2 個あると想定し，各母集団

πi (i = 1, 2) は平均に p 次のベクトル µi ，共分散行列に p 次の正定値対称行列 Σi (> O) をもつと仮定する．ここで，高次元データに対して Σ1 = Σ2 を想定することは現実的ではないので，共分散行列の共通性は仮定しない．ただし，lim inf p→∞ {tr(Σ1 )/tr(Σ2 )} > 0，lim supp→∞ {tr(Σ1 )/tr(Σ2 )} < ∞を仮定する．各母集団 πi から，ni (≥ 2) 個のトレーニングデータ xi1 , ..., xini を無作為に抽出する．
判別対象のデータを x0 (∈ π1 もしくは ∈ π2 ) とし，x0 ∈ πi を誤判別する確率を e(i) と表記する．

∆ = ||µ1 − µ2 ||2 とおく．

高次元の 2 群における判別分析について，Aoshima and Yata (2015) は，∆ → ∞, p → ∞ なる

非スパース性と適当な正則条件のもと，一般的な高次元判別方式に

e(i) → 0, p → ∞, i = 1, 2

(1)なる一致性が得られることを証明した．
一方で，サポートベクターマシン (SVM) は高次元データ解析において疎な解が得られ，汎化性能が良いことも知られているが，Nakayama et al. (2016) は SVM の漸近的性質を高次元小標本の枠組みで研究した．本講演では，トレーニングデータが線形分離可能であることに着目し，線形 SVM(LSVM) の高次元小標本における漸近的性質を導出し，一致性を与えるための正則条件も導出した．
LSVM はある正則条件のもとで (1) を示すが，高次元小標本においてはバイアス項によって不一致性を起こす恐れがある．LSVM はある正則条件のもとで (1) を示すが，高次元小標本においてはバ
イアス項によって不一致性を起こす恐れがある．そこで，そのバイアス項を補正したバイアス補正

LSVM (BC-LSVM) を提案した．さらに，バイアス補正非線形 SVM も提案し，その判別性能を理論的かつ数値的に検証した．

2 高次元データにおける LSVM の漸近的性質とバイアス補正
線形 SVM (LSVM) の判別関数を y(x0 ) とする．いま，p → ∞ で次を仮定する．

(A-i)(A-ii)

Var(||xik − µi ||2 ) = O{tr(Σ2i )}，i = 1, 2;
maxi=1,2 tr(Σ2i )= o(1).
∆21∆∗ = ∆ + tr(Σ1 )/n1 + tr(Σ2 )/n2 とおく．ここで，κ = tr(Σ1 )/n1 − tr(Σ2 )/n2 とし，次を仮定する．

(A-iii)lim sup
p→∞

|κ|
< 1.
∆

定理 1 (Nakayama et al., 2016). (A-i)-(A-iii) を仮定する．p → ∞ のとき判別関数 y(x0 ) について (1) が成り立つ．

(A-iii) は LSVM のバイアス項に関する仮定であることに注意する．しかしながら，通常，高次元の枠組みで n1 と n2 ，もしくは tr(Σ1 ) と tr(Σ2 ) が不均等の場合，(A-iii) は仮定できない．もし，(A-iii) が仮定できない場合は以下の不一致性が成り立つ．
系 1 (Nakayama et al., 2016). (A-i) と (A-ii) を仮定する．p → ∞ のとき判別関数 y(x0 ) について
次が成り立つ．

κ> 1; and∆κe(1) → 0 and e(2) → 1 as p → ∞ if lim sup
< −1.
p→∞ ∆
e(1) → 1 and e(2) → 0 as p → ∞ if lim inf
p→∞

系 1 からもし (A-iii) が成り立たない場合は，もはや LSVM は高次元データ解析において用いるべきではない．そこで，バイアス補正 LSVM (BC-LSVM) を次のように定義する．

yBC (x0 ) = y(x0 ) −

κ̂
.
ˆ
∆∗

ˆ ∗ = ∥x1n −x2n ∥2 ，κ̂ = tr(S 1n )/n1 −tr(S 2n )/n2 である．ただし，xin =
ここで，∆
1
2
1
2
i
と S ini =!nij=1 (xij − xini )(xij − xini )T/(ni − 1)．

!nij=1 xij /ni

定理 2 (Nakayama et al., 2016). (A-i) と (A-ii) を仮定する．p → ∞ のとき判別関数 yBC (x0 ) について (1) が成り立つ．

BC-LSVM は (A-iii) を仮定せずに一致性を示すことができる．
さらに，バイアス補正非線形 SVM も提案し，その判別性能を数値実験と実データ解析を用いて検
証した．

参考文献
[1] Aoshima, M. and Yata, K. (2015). High-dimensional quadratic classifiers in non-sparsesettings. arXiv:1503.04549.
[2] Nakayama, Y., Yata, K. and Aoshima, M. (2016). Support vector machine and its bias
correction in high-dimension, low-sample-size settings. Revised in J. Stat. Plan. Infer.
2Generalized Boltzmann machines
and activation functions
Shinto Eguchi
Institute of Statistical Mathematics, Japan

1

Introduction

Recently the neural network has been revived and attracted lots of interests toward a direction to deep
multilayer networks in machine learning, see LeCun et al. (2015). Boltzmann machines are highly exploredand focused on the general architecture in a paradigm of deep learning. In particular, restricted Boltzmann
machines present a simple understanding for the conditional distributions with logistic sigmoidal functions.
On the other hand, several activation functions including the softplus, leaky rectified linear and max functions rather than the logistic sigmoidal function contribute to eﬃcient learning for network parameters. In
principle, Boltzmann machine is based on Boltzmann distributions which is extended to an exponential
family, cf. Welling (2004) for exponential family harmoniums.This is closely connected with the classicalstatistics that is formalized by the notion of suﬃciency, likelihood, invariance, unbiasedness, eﬃciency andso forth under the assumption of the exponential family.
We present a general framework for Boltzmann machines associated with an activation function, called
the generalized Boltzmann machine (GBM). For this a generalized mean rather than the linear form as in
the arithmetic mean is introduced employing the activation function. We will discuss to combine energyfunctions Eℓ for ℓ, 1 ≤ ℓ ≤ L by the generalized mean#
1 −1 ! 1 "φ
φ(τ Eℓ )
τ
L
L

(1)

ℓ=1

via the activation function φ with the inverse temperature τ . We note that if φ is an identity function
on R, then (1) is nothing but the arithmetic mean; if φ = exp, then (1) equals the log-sum-exponentialfunction. We show that a generator function defines a generalized exponential family. As a result, we givea probabilistic understanding for the GBM including the maxout neural network (Goodfellow et al. 2013) if
we view the generator function as an activation function in a neural network.

2

Generalized Boltzmann machine

We generalize both the Boltzmann distribuiona nd the log-likelihood hunction as follows. Let φ be a strictlyincreasing on R and E an energy function of an input vector x. Then the generalized Boltzmann distribution
is given by%
1 $f (φ) (x, θ) = φ − τ E(x, θ)(2)zθ&with a scale parameter τ and zθ = x φ(−τ E(x, θ)). On the other hand, the generalized loss function for agiven data set D is"1$
L(φ) (θ, D) = −φ−1 f (φ) (x, θ)) + nΨτ (θ),
(3)
τ
D

%
's
& $
where Ψτ (θ) = τ1 x Φ φ−1 (f (φ) (x, θ)) with Φ(s) being −∞ φ(t)dt. We see that the generalization is just to
replace the pair (exp, log) with (φ, φ−1 ), or if φ = exp, then f (φ) (x, θ) is reduced to a Boltzman distribution
and L(φ) (θ, D) equals the minus log-likelihood function up to a constant.
We next consider an energy function of a visible variable x and hidden variable h such that E(x, h, θ) =⊤−b x − c⊤ h − x⊤ W h, where θ = (b, c, W ). In fact E(x, h, θ) is a restricted version of energy function thatis connected only between the visible and hidden components. The generalized Boltzmann distribution is
given by f (φ) (x, h, θ) = φ(−τ E(x, h, θ))/zθ , which has the marginal distribution for the visible variable,f (φ) (x, θ) =1 "φ(−τ E(x, h, θ)).
zθ

(4)

h

Basically we can apply the maximum likelihood for this model (4), however we adopt the generalized estimation method with the loss function!
#" 1
L(φ) (θ, D) = −φ−1 f (φ) (x, θ) + nΨτ (θ).
(5)
τ
x∈D

The gradient is given by a weighted sum of estimating functions as follows.
Theorem 1 Let L(φ) (θ, D) be the objective function defined in (5). Then, the gradient is written by( ∂ φ(−τ E(x, h, θ)) )""∂ (φ)
L (θ, D) =w(φ) (x, h, θ)
∂θ
∂θ
zθ
x∈D h
*"
( ∂ φ(−τ E(x, h, θ)) )+−nEf (φ) (x,θ)w(φ) (x, h, θ)
∂θ
zθ

(6)

h

where
w(φ) (x, h, θ) =φ′ (φ−1 (

&

1
.
φ(−τ
E(x, h, θ))/zθ ))h(7)

We can consider a deep architecture of restricted Boltzmann machines connecting among the input vector
x and hidden vectors h1 , · · · , hL with energy functions.

3

Discussion

The recent developments in deep learning will open a new paradigm for data sciences, however, the understanding to allow uncertainty is not fully elucidated. We discuss to extend Boltzmann machine to permit
various activation functions in a probabilistic argument. This implies to escape from the standard idea by the
exponential family and log-likelihood function. However any extension satisfy a kind of minimax principle.
In eﬀect it is supported if we consider the generalized entropy and divergence. It is necessary to proceedfurther investigations in maximum entropy principle.

References[1] Goodfellow, D. Warde-Farley, M. Mirza, A. Courville & Y. Bengio (2013). ʠ Maxout Networks ʡ(http://arxiv.org/abs/1302.4389). 30th International Conference on Machine Learning.
[2] LeCun, Yann, Yoshua Bengio, and Geoﬀrey Hinton. ”Deep learning.” Nature 521.7553 (2015): 436-444.
[3] Welling, Max, Michal Rosen-Zvi, and Geoﬀrey E. Hinton. ”Exponential family harmoniums with anapplication to information retrieval.” Advances in neural information processing systems. 2004

パラメータ転移学習の理論解析
熊谷亘
神奈川大学 工学部
本研究では, パラメータ転移アプローチを用いた転移学習アルゴリズムについて考察する．特に, パラメトリックな特徴写像に対し，局所安定性とパラメータ転移学習性という新たな概念を導入することで，パラメータ転移アルゴリズムの学習バウンドを導出できることを示す．
従来の機械学習では, データは単一の分布から独立同一に発生すると仮定されている．しかし, この仮定は実際の応用では必ずしも成立しない．そのため, 異なる分布から発生したサンプ
ルを扱うことができる方法を発展させることが重要であると考えられる．このとき，転移学習
はこれらの状況に対処するための一般的な方法を提供する．転移学習では典型的に, 目的のタ
スクに関連する少数のサンプルに加えて, 他のドメインから発生した豊富なサンプルが利用可
能であることが想定されている．そのとき, 転移学習の目的は，他のドメインのデータから有用な知識を抽出し，それを用いて目的のタスクに対するアルゴリズムの性能を向上させることである．転移される知識の種類に応じて, 転移学習の問題を解決するためのアプローチは, イン
スタンス転移, 特徴転移, パラメータ転移などに分類することができる，本研究では, ある種のパラメトリックモデルが想定され，転移された知識はパラメータにエンコードされるようなパラメータ転移アプローチについて考察する．本研究の目的は, パラメータ転移アプローチに基づくアルゴリズムに対して，理論的解析を行うことである．
以下では，パラメータ転移アプローチについて問題設定を述べる．はじめに，幾つかの記法を簡単に導入する．まず，X と Y はそれぞれ，サンプル空間とラベル空間であるとする．
ラベル付きサンプルの空間 Z := X × Y と，その上の結合分布 P (x, y) の組みを領域と呼ぶ．
そのとき，サンプル空間 X と X 上の周辺分布 P (x) の組をドメインと呼び，また，ラベル
集合 Y と 条件付き分布 P (y|x) の組みをタスクと呼ぶ．さらに， H = {h : X → Y} を
仮説空間とし，ℓ : Y × Y → R≥0 を損失関数とする．そのとき, 期待リスクと経験リスクを! n (h) := 1 "n ℓ(yj , h(xj )) によって定める．転移学習の設
R(h) := E(x,y)∼P [ℓ(y, h(x))] と R
j=1n
定において，目標ドメインと呼ばれる，興味のあるドメインから生成されるサンプルに加え，
元ドメインと呼ばれる，別のドメインから生成されるサンプルが利用可能であると仮定される．
本研究では，目標ドメインと元ドメインを区別するために，PT や RS のように，T もしくは
S で表される添字を付すこととする．
以下，目標ドメインでは YT ⊂ R と仮定し, 目標ドメインのパラメトリックな特徴写像 ψθ :1

XT → Rm を用いて，仮説 hT ,θ,w : XT → YT が次のように表されると仮定する:
hT ,θ,w (x) := ⟨w, ψθ (x)⟩.

(1)

ここで，パラメータは θ ∈ Θ と w ∈ WT とし，Θ はノルム ∥ · ∥ が付随したノルム空間の部! T (hT ,θ,w ) を分集合，WT は Rm の部分集合とする．以下では, 単純に RT (hT ,θ,w ) および R!
RT (θ, w) および RT (θ, w) のように記述する. 元ドメインには, サンプル分布 PS,θ,w や仮説
hS,θ,w などのパラメトリックモデルが存在するとし，パラメータ空間の一部 Θ は元ドメインと目標ドメインと共有されていると仮定する．そのとき, θS∗ ∈ Θ と wS∗ ∈ WS は元領域において何らかの指標に関して有効なパラメータであるする．しかしながら, 本研究では θS∗ ∈ Θ とwS∗ ∈ WS に対して明示的な仮定は課さない．
次に，本研究で扱うパラメータ転移アルゴリズムについて説明する．元ドメインと目標ドメ
インでそれぞれ N 個と n 個のサンプルを使用できるとする．パラメータ転移アルゴリズムは，
まず N 個のサンプルを使用して，θS∗ の推定値 θ!N ∈ Θ を出力する．次にアルゴリズムは目標
ドメインのパラメータwT∗:= argmin RT (θS∗ , w)w∈WT

に対し，n 個のサンプルを用いて推定値!T ,n (θ!N , w) + ρr(w)! N,n := argmin Rw
w∈WTを出力する．ここで r(w) は ∥ · ∥2 に関して 1-強凸な正則化項とし，ρ は正の実数とする. 元ドメインが何らかの意味で目標ドメインに関係している場合，元ドメインでの有効なパラメータθS∗ も目標タスクにとって有用であると期待される．RT (θS∗ , wT∗ ) を予測性能の基準値として採用し，学習バウンドを導き出す．局所安定性と転移学習可能性の他に，幾つかの技術的仮定をおくことで，以下の学習バウンドが成り立つ．
定理 1 (学習バウンド). パラメトリック特徴写像 ψθ が局所安定であるとする．また，元ドメ
インで学習された θS∗ ∈ Θ の推定量 θ!N は確率 1 − δ̄ でパラメータ転移学習可能性を満たすとする．ここで， ρ を適切に設定するとき, 以下の不等式が確率 1 − (δ + 2δ̄) で成り立つ :
&%
%
%
%
#
$
1
1
%!
%
%
∗
∗
∗%
!
! N,n − RT (θS , wT ) ≤ C1 √ + C2 %θN − θS % + C3 n 4 %θ!N − θS∗ %. (2)RT θ N , wn
ここで C1 , C2 , C3 は定数である．

サンプル数 N を用いて，推定誤差 ∥θ!N − θS∗ ∥ を評価することができるならば, 定理 1により,どの項が支配的であるかがわかる．特に, 目標ドメインのサンプルと比較して，元ドメインのサンプルがどの程度あれば十分かを評価することができる．

2

大規模な線形予測器のための非同期特徴抽出スキーム

東京大学

情報理工学系研究科

松島 慎

平成 29 年 1 月 27 日
本発表では w ∈ Rp をパラメータとする以下の関数を最小化することを考える：

P (w) ! ∥w∥1 + C

n!

ℓ(⟨w, φ(xi )⟩ , yi ).

(1)

i=1

ここで ℓ : R × Y → R は微分可能な凸関数とし、xi ∈ X と yi ∈ Y は学習に用いる入力と出力とする。また、φ : X → Rp はデータの特徴を記述する関数である。この最小化問題は機械学習の分野で L1 正則化項付きの経験誤差最小化問題と呼ばれる一般的な枠組みで現れ、分類問題や回帰問題だけでなく、構造予測など複雑な出力を入力から予測する場合にも多く用いられる [4]。
一般にサンプル数 n が増大するにつれ、φ として複雑かつ高次元な関数を用いても過学習が起こらず、結果として性能の高い予測器を得ることが可能になる。近年ニューラルネットを用いて表現されるような複雑な予測器を用いて、性能が高い予測が可能になったのも、使用可能なサンプル数が増大したことが一つの要因になっている。このことからも線形予測器においても高い予測器を得るためには、φ として複雑かつ高次元な関数を用いる必要があると考えられる。
標準的に最適化問題 (1) を解く場合、Φ ! {φj (xi )}1≤i≤n,1≤j≤p の計算と実際の求解の段階は分離されている。そのため非常に大きな数の特徴を考えることで Φ のデータ量が計算機のメモリ容量を超す場合に、最適化を効率的に行うことができない問題がある。すなわち、典型的には勾配を求める場合などに Φ の値を読み込む際に、低次の記憶領域にアクセスしなければならず、実際の計算時間が大きくなってしまう。

L1 正則化を用いる上記の問題の場合、最適化問題に必要ない特徴が多く存在し、最適解として疎な解を得られることが期待される。すなわち、予測器が予測に必要な特徴は特徴関数のうち一部の関数のみであり、L1 正則化の影響でこれらに対応するパラメータが零化されることが期待される。そのような場合、結果的に零化される特徴を考えない問題を解くことで実質的に同じ予測器を得ることができる。実際に不必要と推測される特徴を無視して経験的に学習の効率を上げる方法が提案されている。このような方法を利用すれば、必要な特徴のみを展開しメモリに載せることで、上記の問題を回避することができると考えられる。
メモリ容量より大きなデータを用いて正則化付き経験誤差最小化問題を解く最初のスキームとして、Yu

et al. のブロック最小化スキーム [5] が挙げられる。ブロック最小化スキームや Matsushima et al. の Dual
Cached Loops[2] は確率的勾配法を用いることによって (1) の最小化問題を解くことが可能である。しかし、確率的勾配法を用いる場合、最適化問題に必要ない特徴を除外するなど、利用することができない。

1

本発表では、Φ がメモリ容量より大きいような場合に最適化問題 (1) を効率的に解くためのスキームを提案する。本スキームでは二種類の異なるアルゴリズムを非同期的に動作させることで特徴の選択的抽出とパラメータの最適化を同時に行う。本スキームは任意の基底関数について学習することができるスキーム
であるが、基底関数が持つ構造を用いて効率的な学習を行うことも可能である。本スキームでは書き込みス
レッドと訓練スレッドと呼ばれる二つのスレッドがパラメータ w、補助パラメータ u と特徴キャッシュΦJを共有し、非同期的に動作する。書き込みスレッドはデータから特徴 φj を読み込み、これを特徴キャッシュに追加する。ただし、抽出された特徴が現在の解を改良することができなければ、特徴キャッシュに追加し訓練スレッドと共有することなく、展開された特徴は破棄される。すなわち特徴を破棄するための条件は−1 < C∇j L(wt ) < 1.

(2)と記述することができる。この値は前述のように補助パラメータ ut を用いて簡単に計算することができる。
もし特徴キャッシュの容量を超過するようであれば、無作為にキャッシュ内の特徴を破棄し、解放された領域に特徴を追加する。
一方で、訓練スレッドは特徴キャッシュ内に格納されている列に対応する座標を無作為に一つ選び、座標降下法を用いてパラメータ更新を行う。選んだ座標に関して (2) かつ wj = 0 が成り立つ場合、対応する列は特徴キャッシュから破棄することで座標降下法の効率を上げることを目指す。この条件は前節で説明した
ヒューリスティクスに比べより多くの列を不必要とみなすため、実際には必要だった列を破棄する機会がより多くなると考えられる。しかし書き込みスレッドが同時に各列の必要性を確認しているため、より積極的に列を除外することが全体の効率を上げると考えられる。
一般に、書き込みスレッドは特徴を展開することに大部分の時間を費やす一方で、訓練スレッドはパラメータの更新に大部分の時間を費やすため、特徴キャッシュへのアクセスにおける衝突やそれに伴うスレッドの待機はほとんど生じないと考えられる。訓練スレッドに関しては非同期的な座標降下法 [1, 3] を用いて複数のスレッドを利用することができる。また書き込みスレッドは自然に同時に複数のスレッドを利用することができる。さらに特徴空間の構造を利用することにより、より効率的な特徴抽出が可能であると考
えられる。

参考文献
[1] J. Liu, S. J. Wright, C. Ré, V. Bittorf, and S. Sridhar. An asynchronous parallel stochastic coordinate
descent algorithm. Journal of Machine Learning Research, 16(285-322):1–5, 2015.
[2] S. Matsushima, S.V.N. Vishwanathan, and A. J. Smola. Linear support vector machines via dual
cached loops. In Proceedings of Knowledge Discovery and Data Mining, pages 177–185, 2012.
[3] Z. Peng, Y. Xu, M. Yan, and W. Yin. Arock: an algorithmic framework for asynchronous parallel
coordinate updates. SIAM Journal on Scientific Computing, 38(5):A2851–A2879, 2016.
[4] I. Rish and G. Grabarnik. Sparse Modeling: Theory, Algorithms, and Applications. CRC Press, Inc.,2014.
[5] H.-F. Yu, C.-J. Hsieh, K.-W. Chang, and C.-J. Lin. Large linear classification when data cannot fit
in memory. In Proceedings of Knowledge Discovery and Data Mining, pages 833–842, 2010.
2

行列式点過程の準モンテカルロ積分への応用
平尾 将剛（愛知県立大学）∗

はじめに

1

行列点過程 (determinantal point process, DPP) は，1970 年代中頃に Macchi [4] によって量子力学的粒子のひとつであるフェルミ粒子をモデル化するために提案された確率点過程である．この点過程は応用上よく用いられる点過程のひとつであるポアソン点過程とは異なり，反発力がある粒子系を記述することができ，近年盛んに研究されている対象である．
本講演では，近年の行列式点過程を用いた数値積分法への応用について，自身の研究 [3]における結果内容を踏まえ報告を行なった．特に球面アンサンブル，ハーモニックアンサン
ブルと呼ばれる球面上の行列式点過程について着目し，後述するソボレフ空間に対する最悪
誤差の評価を与えた．またその結果，球面アンサンブルは漸近的に最適な収束レートを達成する点列を与えること (準モンテカルロデザイン系列を漸近的に与えること)，及びハーモ
ニックアンサンブルは通常のモンテカルロ法より速い収束レートを達成する点列を与えることを報告した．これらの主結果については次節で必要事項について準備したのちに，3 節においてその詳細を述べる．

準備2

主結果を紹介するため，ここでは幾つかの準備を与える．通常，準モンテカルロ法では，積分領域が d 次元立方体 [0, 1]d のものを考えるが，ここでは d 次元球面 Sd のものを考える．
Brauchart et al. [2] により，近年，球面上の積分に対して準モンテカルロデザイン系列（QMC
design sequence）の概念が導入された．これは球面上の積分に関して，より高速な誤差の収束を実現する球面上の点列を構成する問題である．[2] では，特に球面上に制限された滑ら
かさ s の Sobolev 空間 Hs (Sd ) において準モンテカルロデザイン系列を定義し，その考察を行っている．
次節で与える主結果の証明において重要なのは，Sobolev 空間 Hs (Sd ) は再生核ヒルベ
ルト空間であることである．Brauchart et al. [2] では，内積を適切に選んだ場合，例えば，
d/2 < s < d/2 + 1 に対して，その再生核はK (s) (x, y) := 2Vd−2s (S2 ) − |x − y|2s−d ,x, y ∈ Sd .

と表されることが知られている．また，s > d/2 + 1 においても似た様な表現が得られることが知られている (Brauchart-Womersly [1] や [2] を参照)．
s dSobolev! 空間 H (S ) 上の準モンテカルロデザイン系列を定義するために，球面上の積分"
I(f ) := Sd f (x) dσd (x) と N 点集合 XN ⊂ Sd での近似 Q[XN ](f ) := N1 x∈XN f (x) とにおける最悪誤差 (worst-case error) を次で定義する．

wce(Q[XN ]; Hs (Sd )) :=supf ∈Hs (Sd )∥f ∥Hs ≤1|Q[XN ](f ) − I(f )| .

(1)

このとき，Hs (Sd ) 上の準モンテカルロデザイン系列は次で定義される．
博士 (情報科学) 愛知県立大学 情報科学部（〒 480-1198 愛知県長久手市茨ケ廻間 1522-3, E-mail:hirao@ist.aichi-pu.ac.jp）∗

Definition 2.1 (Hs (Sd ) 上の準モンテカルロデザイン系列, [2]). s > d/2 とし，{XN } を Sd上の点集合の増大列とする．このとき，N に依存しない定数 c(s, d) > 0 が存在し，
wce(Q[XN ]; Hs (Sd )) ≤c(s, d)
.
N s/d

(2)

が成り立つとき，増大列 (XN ) を Hs (Sd ) 上の準モンテカルロデザイン系列であると言う．
特に上の再生核を用いれば，d/2 < s < d/2 + 1 における最悪誤差は，
1 #{wce(Q[XN ]; Hs (Sd ))}2 = Vd−2s (Sd ) − 2|xj − xi |2s−dNi̸=jで与えられることは重要である．主定理ではこの最悪誤差の評価を行なう．

3

主結果

ここで主定理について述べる．我々は前述した 2 種類の球面上の行列式点過程を用いた際の最悪誤差の評価を与える．最初に S2 上の球面アンサンブルを用いた場合である．

Theorem 3.1 ([3]). 1 < s < 2 とする．このとき，N 点球面アンサンブル XN に対して，$%
E wce(Q[XN ]; Hs (S2 ))2 = 22s−2 B(s, N )(3)が成り立つ．ここで B(s, N ) はベータ関数である．

ここでスターリングの公式を用いると，固定した s に対して，B(s, N ) ∼ Γ(s)N −s (N →∞) であることが分かる．したがって，十分大きな N に対して (3) は平均的に (2) を満たしていることが分かる．すなわち，驚くべきことに球面アンサンブルは漸近的に Hs (S2 ) に対する準モンテカルロデザイン系列を与えるのである．
さらに一般次元の球面 Sd に関しては，ハーモニックアンサンブルを用いると次の評価式を得ることができる．この定理は通常のモンテカルロ法より，速いオーダーで最悪誤差が収束していることを述べている．

Theorem 3.2 ([3]). XN を N = dim(PL (Sd )) 点のハーモニックアンサンブルであるとする．このとき，d/2 + 1/2 < s < d/2 + 1 に対して，N とは独立な定数 C(s, d) > 0 が存在し，
E[{wce(Q[XN ]; Hs (Sd ))}2 ] ≤C(s, d)N 1+1/(2d)を満たす．
謝辞 本研究は，科学研究費補助金（若手研究 (B)，課題番号 16K17645）の助成を受けている．

References[1] J.S. Brauchart, R.S. Womersley. Numerical integration over the unit sphere, L2 discrepancy and sum of distances. In preparation.
[2] J.S. Brauchart, E.B. Saﬀ, I.H. Sloan, R.S. Womersley. QMC Designs: optimal order
quasi-Monte Carlo integration schemes on the sphere. Math. Comp., 83(290): 2821–
2851, 2014.
[3] M. Hirao. QMC designs and determinantal point processes. Submitted to MCQMC2016 conference proceedings.
[4] O. Macchi. The coincidence approach to stochastic point processes, Adv. Appl. Prob.,
7: 83–122, 1975.

Cell Regression and Reference Priors
Jinfang Wang and Shigetoshi Hosaka
Chiba University and Hosaka Clinic of Internal Medicine

1. Introduction Many data collected by survey agencies, such as government offices, are published in table forms. A large portion of these tables show
cell frequencies for categorized continuous variables. In this paper, we shallrefer to such coarse tables as reference tables. Cell regression methods proposedin this paper only explore the information on these cell frequencies and theintervals defining the cells. We use parametric models to predict the cell probabilities, and estimate the regression parameters by minimizing the KullbackLeibler divergence from the reference distribution to the predictive distribution. Bayesian extensions are also considered. We propose an approximate
Markov chain Monte Carlo algorithm to compute the posterior distribution of
the parameters. The posterior distribution so obtained may be used as priordistribution in a second stage analysis based on more detailed data. For thisreason, we call the posterior distribution the reference prior. Indeed, the primary
motivation of cell regression is to mining prior information embedded in often
large but coarse tables available with low costs.
2. Predictive Distributions We assume that both xi and yi are continuous,and are interval-censored. The exact values of xi and yi are not known. We
only know the cell membership for each pair ( xi , yi ), for i = 1, · · · , R.
Definition 1 (Predictive Cell Distributions). Assume that the following cell regression model holds.
yi |( xi , θ) ∼ f (yi | xi , θ)xi ∼ f ( xi )

(1)
(2)

where both f (yi | xi , θ) and f ( xi ) are known up to an unknown parameter θ. A predictive cell distribution pij (θ) is defined by#! c j " ! dipij (θ) =f ( x, y|θ) dy dx(θ ∈ Θ) .
(3)
c j −1d i −1Let Rij be cell frequencies, rij = Rij /R and R = ∑ Rij . Let P be the predictive distribution defined by (3). The goodness of P is measured through the
Kullback-Leibler divergence from the reference distribution R to P ,KL(R||P ) = ∑ rij log
i,j

1

rij
pij

(4)

We define the minimum contrast estimator $θ to be the solution to the following
constrained minimization problem,minθ

KL(R||P )subject to∑ pij (θ) = 1

(5)

i,j

3. Reference Priors Bayesian random effects models are often used in deriving the posterior predictive distribution for a particular subject based on
repeated measurements for each subject. The cell regression analyses comeinto play when we want to replace the assumption on the prior distribution of
the fixed parameter. We propose to use the posterior distribution of the corresponding parameters derived from a Bayesian cell regression analyses as theprior distribution for the fixed parameter. We consider the following Bayesiancell regression model.
yi |( xi , θ) ∼ f (yi | xi , θ)(6)

θ ∼ π (θ)

(8)

(7)

xi ∼ f ( xi )We define the Bayesian predictive cell distribution pijb (θ) by#! c j " ! dibpij (θ) =f ( x, y) dy dx(θ ∈ Θ)It is easily seen thatc j −1(9)d i −1pijb (θ) = π (θ) pij (θ)(10)

Denote by P b the Bayesian predictive distribution. Letb$θ = arg min KL(R||P b ) = arg min KL(R||P ) − log π (θ)
θ∈Θ

θ∈Θ

(11)

The following algorithm computes the posterior distribution for θ.
b(i) Compute $θ of (11).
(ii) For the (i, j) cell, generate Rij independent pairs of ( x, y) according tox ∼ f (x)(12)

by |( x, θ) ∼ f (y| x, $
θ)

(13)

And do this for all cells to generate R independent samples( x1 , y1 ), · · · , ( x R , y R )(14)

yi |( xi , θ) ∼ f (yi | xi , θ)(15)

θ ∼ π (θ)(17)(iii) Apply a typical Markov chain Monte Carlo method to compute the posterior distribution of θ based on the approximate sample (14) and theBayesian cell regression modelxi ∼ f ( xi )

2

(16)

Euclidean Design Theory
澤正憲
神戸大学大学院システム情報学研究科
sawa@people.kobe-u.ac.jp
関数の積分値を，定義域の有限個の点 (ノード) での関数値の重み付き平均で近似する公式を，cubature 公式 (cubature formula) という．特に，次数 t以下の任意の多項式 f について，積分値と f の重み付き平均値の誤差が 0 になるとき，次数 t の公式 (cubature of degree t) という．

Simpson や Gauss の名を冠する公式があるように，少ないノードからなり，かつ次数の大きな公式の存在および構成問題は，直交多項式の零点の解析と並行して，数値解析，特殊関数論などの分野を中心に古くから調べられてきた．実は，そのような cubature 公式は• 実・複素幾何における球面デザイン (spherical design)，• 微分幾何学における大対蹠集合 (great antipodal set)，• 離散幾何学における距離集合 (distance set)，• 関数解析学におけるバナッハ空間のノルム不変線形作用素 (isometricembedding)，• 高次形式論におけるヒルベルト恒等式 (Hilbert identity)，• 整数論のテータ級数に関する Lehmer 予想 (Lehmer conjecture)，• 量子情報理論における mutually unbiased basis (MUB)，• 統計的実験計画法における D 最適計画 (D-optimal design)，などの異分野の対象と深く結び付いているのだが，その多くは個々の分野で独自に調べられてきた経緯があり，最近でも，双方にとって有益な研究成果が互いに認識されていないというケースも少なくない．
そこで，本研究では，上述の諸研究対象を統一的に扱うべく，被積分関数のクラスを一般の関数空間に押し上げて，測度空間上の cubature 公式 (cubature

on measure space) の概念を導入した．そして，以下の主成果を得た：1

(1) 多項式型 cubature 公式のノードの個数に関する Stroud 不等式 (Stroudbound) の一般化．特に BIB デザインのブロック数に関する Fisher 不等式の解析的別証明の提示．
(2) (1) の不等式においてタイト (tight) な公式の再生核による特徴付け．
(3) 多項式型 cubature 公式の漸近存在定理として知られている Tchakaloﬀの定理 (Tchakaloﬀ Theorem) の一般化．特に位相幾何学的な諸条件下での Chebyshev 型の cubature 公式，すなわち「重み付き平均」が「算術平均」であるような公式の漸近存在証明．

(4) 多項式型 cubature 公式に対する Sobolev の不変式論の一般化．
(5) 多項式型 cubature 公式の逐次的構成法の一般化．
このように，「測度空間上の cubature 公式論」では，多項式型の cubature
公式に関する諸概念・諸定理を一般化して，異分野の研究対象を網羅的に扱うための理論的な枠組みを提供する．また，統計的には，BIB デザインの Fisher
不等式の解析的別証明 (主結果 (1))，直交配列の「直交性」の解析的な解釈，
そして D 最適計画の構成の提示など，主に実験計画法への応用に期待感がある．
本研究は拙論文 [2] に基づいている．なお，多項式型の cubature 公式の理論と応用については，Dunkl-Xu [1], Sobolev-Vaskevich [3], Stroud [4] などに詳しく書かれているので，そちらを参照されるとよい．そのうち，測度空間上の cubature 公式論の統計学者向けの書「Euclidean Design Theory」がSpringer Briefs から出版される予定なので，そちらも手にとっていただけると嬉しく思います．

参考文献[1] C.F. Dunkl, Y. Xu. Orthogonal Polynomials of Several Variables.
Encyclopedia of Mathematics and its Applications, 81. Cambridge University Press, 2001.
[2] M. Sawa. Cubature 公式の理論．数学 62 (2016), 24-53.
[3] S.L. Sobolev, V.L. Vaskevich. The Theory of Cubature Formulas.
Mathematics and its Applications, 415. Kluwer Academic Publishers
Group, Dordrecht, 1997.
[4] A.H.

Stroud. Approximate Calculation of Multiple Integrals.

Prentice-Hall, N.J., 1971.

EXACT VC DIMENSION OF ELLIPSOIDS AND CONSISTENCYOF MAXIMUM LOG-LIKELIHOOD ESTIMATOR FOR
MULTIVARIATE GAUSSIAN MIXTURES
YOHJI AKAMA

1. Introduction
The log-likelihood function Ln (θ) to estimate the parameter θ of a d-dimensional(d > 1) K-component Gaussian mixture is unbounded, where n is the size of a sample. This is overcome in Chen and Tan [4], by adding a term pn (θ) to Ln (θ). This
penalizes small variances or ratios of variances to Ln (θ). Chen and Tan attempted
to prove the consistency of their penalized maximum likelihood estimator (MLE),
by reducing the argument to a univariate argument combined with a Bernsteininequality and Borel-Cantelli Lemma. But, at some point, they did not uniformly
handle the unbounded subsets of Rd , although they should. See [3], for detail. To fix
the flaw in Chen and Tan’s consistency proof, Alexandrovich slightly modified Chen
and Tan’s suﬃcient condition for the penalized MLE to be consistent, and used (1)a uniform law of iterated logarithm for classes having finite VC dimensions and (2)that the set of the d-dimensional ellipsoids has finite VC dimension. VC dimensionof a class is a combinatorial measure of the complexity of the class. It is diﬃcult
to give the exact value, in general. In this note, we give the exact VC-dimensionof the d-dimensional ellipsoids (Akama-Irie [1]). Then we review Alexandrovich’s
approach to consistency proof [3] of Chen and Tan’s penalized MLE to Gaussian
mixture.
For sets X ⊂ Rd and Y ⊂ X, we say that a set B ⊂ Rd cuts Y out of X if
Y = X ∩ B. A class C of subsets of Rd is said to shatter a set X ⊂ Rd , if for any
Y ⊂ X there exists B ∈ C such that B cuts Y out of X.
The Vapnik-Chervonenkis dimension (VC dimension for short) of a class C is
VCdim(C) := sup{ ♯S | S ⊂ Rd , C shatters S }.

By a d-dimensional ellipsoid, we mean a region in Rd defined as {x ∈ Rd |Q(x − µ) < 1}, where µ ∈ Rd and Q is a positive definite quadratic form. Our main
result is:
Theorem 1. The class of d-dimensional ellipsoids has VC dimension (d2 + 3d)/2.
2. Alexandrovich’s consistency proof of Chen and Tan’s penalized
MLE
In order to uniformly handle unbounded subsets of Rd missed in Chen and Tan’sconsistency proof, Alexandrovich slightly modified the condition C3 , so that he can
employ uniform law of iterated logarithm developed in empirical process theory [5].
Condition 1.
C̃3 : p̃n (Σ) ≤ a(n) log |Σ|for |Σ| ≤ cn−2d . Here c is a constant and a(n) = o(n).
12YOHJI AKAMA!K
Here n is the size of a sample, pn (θ) =j=1 p̃n (Σj ) is the penalty function
for a K-component Gaussian mixture parameter θ = (µ1 , Σ1 , p1 , . . . , µK , ΣK , pK ),!Kµj ∈ Rd , a positive definite Σj ∈ Rd×d , and j=1 pj = 1.
From Chen and Tan’s argument [4], Alexandrovich observed that the following
condition implies the consistency of Chen and Tan’s penalized MLE, if the penalty
function satisfies C̃3 instead of C3 . Condition 2 gives an upper bound of the numberof the sample random variables lying in a ellipsoid, asymptotic in the size n of thesample and the “volume” of the ellipsoid.
Condition 2. Yn (n ∈ N) is an i.i.d. process such thatn"1{(Yi −µ)⊤ Σ−1 (Yi −µ)≤(log |Σ|)2 } ≤ a(n) + b(n, |Σ|)i=1

b(n,s)
where a(n) = o(n), b(n, s) = O(n), lims→0+ √
= 0.
log s

Proposition 1 (Alexander [2]). Suppose that C is a class of Borel subsets of Rd
and C has a finite VC dimension. Let Yn (n ∈ N) be a d-dimensional i.i.d. process.
Then a.s.,!n#| i=1 1C (Yi ) − nPY1 (C)|√lim sup sup= sup PY1 (C)(1 − PY1 (C)).
2n log log nn→∞ C∈CC∈C
By Theorem 1 and Proposition 1 with C being the set of d-dimensional ellipsoids,
Alexandrovich derived:

Corollary 1. Let (Yn )n∈N be a d-dimensional i.i.d. process with a bounded Lebesgue
density f , M := supy f (y). Then a.s. there exists a positive integer N such thatn"i=1

1{(Yi −µ)⊤ Σ−1 (Yi −µ)≤(log |Σ|)2 } ≤3#n log log n + nM vd |Σ|1/2 (log |Σ|)d4

for every µ ∈ Rd , every positive definite matrix Σ ∈ Rd×d , and every n ≥ N . Herevd is the volume of the d-dimensional unit ball.
Then this satisfies Condition 2. Alexandrovich’s observation explained above
implies the consistency of Chen and Tan’s penalized MLE under the modified condition C̃3 .
References[1] Akama, Y. and Irie, K., VC dimension of ellipsoids, arXiv:1109.4347 [math.CO]., (2011)[2] Alexander, K. S., Probability inequalities for empirical processes and a law of the iterated
logarithm, Annals of Probability 12, 1041―1067 (1984).
[3] Alexandrovich, G., A Note on the Article ’Inference for multivariate normal mixtures’ by J.
Chen and X . Tan., J. Multivariate Analysis 129, 245 ― 248, 2014.
[4] Chen, J. and Tan, X., Inference for multivariate normal mixtures, Journal of Multivariate
Analysis, 100 (2009), 1367–1383.
[5] Dudley, R.M., Uniform central limit theorems, Cambridge Studies in Advanced Mathematics,
vol. 63, Cambridge University Press, 1999. MR MR1720712 (2000k:60040)
Mathematical Institute, Tohoku University, Sendai, 980-8578, Japan
E-mail address: akama@m.tohoku.ac.jp

入れ子型混合モデルに基づく cancer outlier profile の推定とがん診断への応用 について
松井 孝太 1 , 大浦 智紀 2 , and 松井 茂之 1
1

名古屋大学大学院医学系研究科 / JST CREST
2
Eli Lilly Japan

1 はじめにがん患者の集団における発がん遺伝子の異質性は, 疾患生理学の理解, リスクグループの特定および患者の治療の最適化に重要な意味を有すると考えられている. Tomlin et al. [4] は, 前立腺がんを対象とした研究で , 従来の 2 標本 t-検定のような全がん患者で共通な発現プロファイルを持つ遺伝子を特定しようとする手法では捉えられないプロファイルが存在することを示した. そのような,
一部のがん患者において一部のがん遺伝子が示す特異的な発現プロファイルを outlier expression
profile (OEP) と呼び, そのようながん遺伝子のことを cancer outlier と呼ぶ. cancer outlier は, 上述した異質性の代表的な例と考えられ, [4] では OEP を持つ遺伝子の特定する方法として cancer outlier
profile analysis (COPA) が開発され, これに基づいて前立腺がん患者のサブタイプが明らかとなった.
Tomlin et al. の研究に触発された COPA の研究が, その後複数のグループによって実施されている. Wu [5] は, 健常者の発現プロファイルを基準として cancer outlier を定義し, これを特定するための新しい検定統計量として ORT 統計量を提案した. Tibshirani and Hastie [3] は, Wu とは異なる
cancer outlier の定義に基づいて outlier-sum statistic を用いた検定による cancer outlier のスクリー
ニングを提案した. Tomlin et al. を含めたこれらの COPA 研究は, それぞれ独自に定義した cancer
outlier に対して行われている. すなわち, 健常者のプロファイルからどの程度乖離していれば outlierと見做すのか, というしきい値が研究ごとに異なる. Lian [1] は, これらの問題設定を統一的に評価するため, すべてのしきい値に関して考察を行い, MOST statistic を用いた outlier 遺伝子の検定手法を開発した.
ここまでに紹介した COPA 研究は, すべて cancer outlier として振舞っている遺伝子を仮説検定に基づいて特定するための手法の研究である. しかし, これらの多重検定を用いた outlier 遺伝子のスクリーニングは, cancer outlier analysis の第一歩に過ぎない. 例えば, がん関連遺伝子のあるグループは, 発がん経路で共調節されたり, また増幅された染色体領域や遺伝子融合に関連する染色体領域に共局在する特性を有する. したがって, よく似た OEP を持つがん関連遺伝子をクラスタリ
ングすることで, このような生物学的にも意味のある遺伝子群の特定に繋がる. さらに, そうして特定された遺伝子群を用いた, cancer outlier に基づくがん患者の識別法 (すなわち, がんの診断法) の開発も期待される.
以上のような背景の下で, 本研究ではパラメトリック入れ子型混合モデルに基づいた遺伝子とサンプル両方向のクラスタリング (biclustering) 手法を提案する (Figure 1). 提案法は, cancer outlier
を陽にモデル化し, さらに表現型 (がん or 健常) の情報も取り入れた教師付きのクラスタリング
手法となっており, EM アルゴリズムを用いてパラメータを推定する. また, 適切な OEP 遺伝子のコンポーネント数を決定するために, 情報量基準に基づくモデル選択を行う. 一度提案法によってFigure 1 のような構造 (outlier 遺伝子の従う分布) が推定されれば, 新たな患者に対して, がんか健常かを確率的に判別するアルゴリズムを構成することができる. 本研究では, outlier 構造の推定か
らがん患者の判別までの一連のプロシージャを提案法として, シミュレーション及び実データによ1SampleCancer
Cancerrelatedgenes

Gene

NormalNonrelated
genes
Figure 1: 本研究で提案する biclustering の結果の概念図. 行方向 (遺伝子方向) に見ると, がん関連
遺伝子 (青の部分) と関連なし遺伝子 (白の部分) に分かれる. 列方向 (患者方向) に見ると, がん患者の群において, 関連遺伝子の発現プロファイルが患者毎に異なる.

る評価実験を行う. 実データとして, Mils et al. [2] で解析された骨髄異形成症候群の表現型データ
及びマイクロアレー遺伝子発現量データを用いる.

References[1] Heng Lian. Most: detecting cancer diﬀerential gene expression. Biostatistics, 9(3):411–418, 2008.
[2] Ken I Mills, Alexander Kohlmann, P Mickey Williams, Lothar Wieczorek, Wei-min Liu, Rachel Li,
Wen Wei, David T Bowen, Helmut Loeﬄer, Jesus M Hernandez, et al. Microarray-based classifiers and prognosis models identify subgroups with distinct clinical outcomes and high risk of aml
transformation of myelodysplastic syndrome. Blood, 114(5):1063–1072, 2009.
[3] Robert Tibshirani and Trevor Hastie. Outlier sums for diﬀerential gene expression analysis. Biostatistics, 8(1):2–8, 2007.
[4] Scott A Tomlins, Daniel R Rhodes, Sven Perner, Saravana M Dhanasekaran, Rohit Mehra, XiaoWei Sun, Sooryanarayana Varambally, Xuhong Cao, Joelle Tchinda, Rainer Kuefer, et al. Recurrent
fusion of tmprss2 and ets transcription factor genes in prostate cancer. science, 310(5748):644–648,
2005.
[5] Baolin Wu. Cancer outlier diﬀerential gene expression detection. Biostatistics, 8(3):566–575, 2007.

2࣮༻తͳՃ଎ۙ઀ޯ഑๏ͷ࣮૷ͱ 2 ஋൑ผϞσϧͷԠ༻*ҏ౻௚‫( ل‬౦‫ژ‬େֶ), ෢ా࿕ࢠ (౷‫਺ܭ‬ཧ‫)ॴڀݚ‬, TOH Kim-Chuan (γϯΨϙʔϧࠃཱେֶ)2017 ೥ 2 ݄ 19 ೔ Պ‫ݚ‬අγϯϙδ΢Ϝʮ౷‫ܭ‬తϞσϦϯάͱ‫ࢉܭ‬ΞϧΰϦζϜͷ਺ཧͱల։ʯ
ͨͩ͠, Lk+1 = Lk ≥ Lf (k = 1, 2, . . .) ͱ͓͘. ࣮͸,F (α) Λ β ͷ·ΘΓͰۙࣅͨؔ͠਺:1 ͸͡Ίʹ
2 ஋൑ผ໰୊͸‫ػ‬ցֶशʹ͓͚Δॏཁͳ໰୊ͷ 1 ͭͰ͋Γ, αϙʔτϕΫλʔϚγϯ (SVM) Λ͸͡Ίͱͯ͠
༷ʑͳ 2 ஋൑ผϞσϧ͕ఏҊ͞Ε͍ͯΔ. ߴ͍൑ผਫ਼౓Λୡ੒͢ΔͨΊʹ͸, σʔλ͝ͱʹ༷ʑͳϞσϧΛղ͖,ద੾ͳϞσϧΛબ୒͢Δඞཁ͕͋Δ. ͦͷͨΊ, ࣗ༝ͳϞσϧબ୒Λߦ͏্Ͱ͸, ‫ݸ‬ʑͷϞσϧʹಛԽͨ͠ղ๏
Ͱ͸ͳ͘, ൚༻త͔ͭߴ଎ͳ࠷దԽख๏Λ։ൃ͢Δ͜ͱ
͕ॏཁͰ͋Δ. ͦ͜ͰຊߘͰ͸ҎԼͷ࠷దԽ໰୊ʹର͢Δղ๏Λߟ͑Δ.
min F (α) := f (α) + g(α).

QL (α; β) = f (β) + ⟨∇f (β), α − β⟩ + g(α) + L2 ∥α − β∥22Λ༻͍Δͱ,#$
proxg,Lk αk − L1 ∇f (αk ) = argminα∈Rd QLk (α; αk )ͱද͢͜ͱ͕Ͱ͖Δ. ͢ͳΘͪ, PG ͸, ؔ਺ F ͷۙࣅؔ਺ QLk (α; αk ) Λ‫܁‬Γฦ͠࠷খԽ͢ΔΞϧΰϦζϜͩͱ‫ͤͳݟ‬Δ. ͜ͷ PG Ͱੜ੒͞ΕΔ఺ྻʹ “‫׳‬ੑ” ͷΑ͏ͳಈ͖ΛՃ͑ͨ΋ͷ͕, Ճ଎ۙ઀ޯ഑๏ (APG) [1] Ͱ
͋Δ. APG Ͱ͸ҎԼͷ༷ʹ఺ྻΛߋ৽͢Δ.

(1)

α∈Rd

#
$
Step 1: αk ← proxg,Lk β k − L1k ∇f (β k ) .
√1+ 1+4t2k
Step 2: tk+1 ←.
(ͨͩ͠ t1 = 1.)
2
tk −1
k+1
k
k
k−1
Step 3: β← α + tk+1 (α − α).

ͨͩ͠, ҎԼͷԾఆΛ͓͘.

• f : Rd → R ͸࿈ଓతඍ෼Մೳͳਅดತؔ਺Ͱ, ͦͷޯ഑ ∇f (·) ͸Ϧϓγοπ࿈ଓͰ͋Δ. ͢ͳΘͪ,͋ΔϦϓγοπఆ਺ Lf > 0 ͕ଘࡏͯ͠, ҎԼͷෆ౳͕ࣜ੒Γཱͭ.

PG ͷऩଋϨʔτ͸ F (αk ) − F (α∗ ) ≤ O(1/k) Ͱ͋Δͷʹରͯ͠, APG ͸ΑΓ଎͍ऩଋϨʔτ F (αk ) −F (α∗ ) ≤ O(1/k 2 ) Λ΋ͭ. FAPG [2] ͸, APG ͱ༷ʑͳ࣮༻্ͷվྑ๏Λ૊Έ߹Θͤͨख๏Ͱ͋Γ, Algorithm 1 ͷΑ͏ʹ఺ྻΛߋ৽͢Δ. FAPG ͷऩଋϨʔτ%&͸ F (αk ) − F (α∗ ) ≤ O (log k/k)2 ͕ͩ, ࣮༻తʹ͸
APG ΑΓ΋ߴ଎Ͱ͋Δ.

∥∇f (α) − ∇f (β)∥2 ≤ Lf ∥α − β∥2 ∀α, β ∈ Rd .
• g : Rd → R ∪ {+∞} ͸ਅดತؔ਺Ͱ, ༗ޮྖҬdom(g) = {α ∈ Rd | g(α) < +∞} ͸ดತͰ͋Δ.
‫ػ‬ցֶश΍৴߸ॲཧͰදΕΔ༷ʑͳϞσϧ͕, f Λଛࣦ߲, g Λਖ਼ଇԽ߲ͱஔ͘͜ͱͰ໰୊ (1) ʹ‫ؼ‬ணͰ͖Δ.
໰୊ (1) ʹର͢Δޮ཰తͳղ๏ͱͯ͠, Ճ଎ۙ઀ޯ഑๏

3 ઢ‫ܗ‬൑ผϞσϧ΁ͷԠ༻(APG) [1] ͕஌ΒΕ͓ͯΓ, APG Λվྑ࣮ͯ͠༻্ͷੑೳΛߴΊΔ‫ڀݚ‬΋ͳ͞Ε͍ͯΔ. Ұ෦ͷվྑ๏ʹରͯ͠͸େҬతऩଋੑ͕ࣔ͞Ε͍ͯͳ͔͕ͬͨ, զʑ [2] ͸ෳ਺ͷվྑ๏Λ૊Έ߹ΘͤΔ͜ͱͰ, ཧ࿦తʹऩଋϨʔτͷอূͷ͋Δ࣮༻తʹߴ଎ͳ APG (FAPG) ΛఏҊ͠
ͨ. ຊߘͰ͸ FAPG ͷద༻ྫΛ঺հ͢Δ.

2 ஋൑ผ໰୊͸, ಛ௃ྔ x ͱϥϕϧ y ͷ m ‫ݸ‬ͷ૊:(xi , yi ) ∈ Rn × {+1, −1}, i ∈ M := {1, . . . , m} Ͱ༩͑ΒΕͨ‫܇‬࿅σʔλΛ΋ͱʹ, ະ஌ͷσʔλ x̂ ͕ͲͪΒͷΫϥεʹଐ͍ͯ͠Δ͔Λ൑ผ͢Δ໰୊Ͱ͋Δ. ͍·, ઢ‫ ਺ؔܗ‬h(x̂) = w ⊤ x̂ ͷූ߸ʹΑͬͯ σʔλ x̂ ͷΫϥεΛ༧ଌ͢Δ͜ͱΛߟ͑Δ. ͜ͷͱ͖ద੾ͳ܎਺ w ͷఆΊํΛ࠷దԽ໰୊ͱͯ͠ఆࣜԽͨ͠΋ͷ͕ઢ‫ܗ‬൑ผϞ

2 ఏҊΞϧΰϦζϜσϧͰ͋Δ. ༷ʑͳઢ‫ܗ‬൑ผϞσϧ͕໰୊ (1) ʹ‫ؼ‬ணؔ͞਺ g ͷۙ઀ࣸ૾͸ҎԼͷ༷ʹఆٛ͞ΕΔ.

!

Ε, FAPG ʹΑͬͯޮ཰తʹղ͘͜ͱ͕Ͱ͖Δ.

"

3.1 ℓ1 ਖ਼ଇԽϞσϧg(w) = ∥w∥1 ͱ͓͘. ͜ͷͱ͖,'m• f (w) = C i=1 log(1 + exp(−yi w⊤ xi )) ͱ͓͘ͱ,໰୊ (1) ͸ ℓ1 ਖ਼ଇԽϩδεςΟοΫճ‫ͳͱؼ‬Δ.
'm• f (w) = C i=1 (max{0, 1 − yi w⊤ xi })2 ͱ͓͘ͱ,໰୊ (1) ͸ ℓ1 ਖ਼ଇԽ ℓ2 -SVM ͱͳΔ.

proxg,L (ᾱ) = argmin g(α) + L2 ∥α − ᾱ∥22 .
α∈Rdۙ઀ޯ഑๏ (PG) ͸, ؔ਺ f ͷޯ഑ͱ, ؔ਺ g ͷۙ઀ࣸ૾Λ༻͍ͯ, ҎԼͷ༷ʹ఺ྻΛߋ৽͢Δ.

#$
αk+1 ← proxg,Lk αk − L1k ∇f (αk ) .

1

U ′ ͱ x(α) ͷఆΊํͱ, 2 ஋൑ผϞσϧͷؔ࿈. ͨͩ͠, M+ = {i ∈ M | yi = +1}, M− = {i ∈ M |yi = −1} Ͱ͋Δ. ·ͨ, αϯϓϧू߹ {xi | i ∈ Mo } ʹରͯ͠, x̄o ͸ฏ‫ۉ‬ϕΫτϧ, Σo ͸‫ڞ‬෼ࢄߦྻΛද͢(o ∈ {+, −}).
ද1

(MPM)
ϑΟογϟʔͷઢ‫ܗ‬൑ผ‫ث‬(FDA)
ν-αϙʔτϕΫλʔϚγϯ(ν-SVM)U ′ = {α ∈ Rn | ∥α∥2 ≤ κ}
%
& %
&1/2
x(α) = x̄+ − x̄− + Σ+ + Σ−
α
'
'
1
′
m
U = {α ∈ R | i∈M+ αi = i∈M− αi = 12 , 0 ≤ α ≤ mνe}''x(α) = i∈M+ αi xi − i∈M− αi xiͨͩ͠, w = x(α)/∥x(α)∥2 Ͱ͋Δ. ද 1 ͸ U ′ ͱ x(·)

Algorithm 1 FAPG Method0ͷఆΊํͱ, 2 ஋൑ผϞσϧͷؔ࿈Λද͍ͯ͠Δ. ͍·,

Input: L0 > 0, ηu , ηd > 1, α , K1 ≥ 2, δ ∈ (0, 1)
Init.: t1 ← 1, t0 ← 0, β 1 ← α0 , α−1 ← α0
Init.: L1 ← L0 , i ← 1, kre ← 0
for k = 1, 2, . . . do!
"
αk ← proxg,Lk β k − L1 ∇f (β k )# Step 1

f (α) = ∥x(α)∥2 , g(α) = δU ′ (α) ͱ͓͚͹, FAPG ʹΑͬͯ (3) Λղ͘͜ͱ͕Ͱ͖Δ. ͍ͣΕͷྫʹ͓͍ͯ΋,ۙ઀ࣸ૾ proxg,L (w) ͸ Breakpoint ୳ࡧ๏ [3] ͳͲΛ༻͍ͯઢ‫ ؒ࣌ܗ‬O(m) Ͱ‫͢ࢉܭ‬Δ͜ͱ͕Ͱ͖Δ.

kwhile F (αk ) > QLk (αk ; β k ) do
Lk ← ηu!
Lk
1+

1+4(Lk /Lk−1 )t2
k−1

4 ਺஋࣮‫ݧ‬# ‘bt’tk ←2
t−1
β k ← αk−1 + k−1(αk−1 − αk−2 )! tk"
αk ← proxg,Lk β k − L1k ∇f (β k )end while# ‘dec’# ‘dec’105# ‘bt’

103

Lk+1 ← Lk√/ηd# ‘dec’1+

1+4(L/L )t2

k+1
k k
tk+1 ←
2
tk −1
k+1
k
k
β
← α + tk+1 (α − αk−1 )Time (sec.)

104

Time (sec.)

2

105
104

102
103
Number of features n104

m=10000, C=1
LIBLINEAR
FAPG

103
102
101
100

104
Number of samples m

105

10 1 110102
103
Number of features n104αϯϓϧ਺ m ͓Αͼ࣍‫ ݩ‬n ͷେ͖͞ͱ‫ؒ࣌ࢉܭ‬ͷؔ܎. (্ஈ: ν-SVM. Լஈ: ℓ1 ਖ਼ଇԽ ℓ2 -SVM.)

FAPG Λద༻͢Δ͜ͱ͕Ͱ͖Δ. ͜ͷ g(w) ͷ΋ͱͰͷۙ઀ࣸ૾͸ιϑτᮢ஋ؔ਺ͱ‫ݺ‬͹Ε, (proxg,L (w))i =sign(wi ) max{0, |wi | − L} ͱ‫͖Ͱࢉܭ‬Δ.
3.2 ౷Ұత൑ผϞσϧℓ2 ਖ਼ଇԽϞσϧΛ‫ؚ‬Ή༷ʑͳઢ‫ܗ‬൑ผϞσϧ͕, ҎԼͷ‫Ͱࣜܗ‬౷Ұతʹ‫ه‬ड़͢Δ͜ͱ͕Ͱ͖Δ.

ν-SVM ͱ ℓ1 ਖ਼ଇԽ ℓ2 -SVM ʹର͠, ਓ޻σʔλΛ༻͍ͯ‫ط‬ଘͷιϑτ΢ΣΞͱ FAPG ͷ‫ؒ࣌ࢉܭ‬Λൺֱ͠
ͨ. ಛʹେ‫ن‬໛ͳσʔλΛ༻͍ͨͱ͖ʹ, FAPG ͷ༏Ґੑ͕‫ݟ‬ΒΕͨ.

ࢀߟจ‫ݙ‬(2)[1] A. Beck and M. Teboulle: A fast iterative shrinkagethresholding algorithm for linear inverse problems.
SIAM Journal on Imaging Sciences, 2(1):183–202,
2009.
[2] N. Ito, A. Takeda, and K.-C. Toh: A unified formulation and fast accelerated gradient method for classification. Journal of Machine Learning Research, toappear.
[3] K.C. Kiwiel: Breakpoint searching algorithms for thecontinuous quadratic knapsack problem. Mathematical Programming, 112:473–491, 2008

ͨͩ͠, φU (w) = maxx∈U {−w ⊤ x} Ͱ͋Γ, U ⊆ Rm

͸ดತू߹Ͱ͋Δ. ू߹ U ͷେ͖͞΍‫ܗ‬Λద੾ʹఆΊΔ͜ͱʹΑΓ, ༷ʑͳ 2 ஋൑ผϞσϧ͕ (2) ʹ‫ؼ‬ணͰ͖Δ.

͍·, ഔհม਺ α ∈ Rd , ઢ‫ ૾ࣸܕ‬x(·) : Rd → Rn , ดತू߹ U ′ ⊆ Rd Λ༻͍ͯ, U = {x(α) | α ∈ U ′ } ͱද͢ͱ, ໰୊ (2) ͸ҎԼͷ࠷খϊϧϜ໰୊ʹ‫ؼ‬ண͞ΕΔ.
α∈U

n=1000, C=1
LIBLINEAR
FAPG

101

ͨͩ͠ C > 0 ͸ύϥϝʔλͰ͋Δ. ͍ͣΕͷ৔߹΋

max min −w⊤ x(α) = − min′ ∥x(α)∥2 ,SeDuMiLIBSVM
FAPG
LIBLINEAR100

10 3 1
10

105

102

ਤ1

101

10
104
Number of samples m

103

100 3
10

102

10 1

105

end forα∈U ′ ∥w∥2 ≤1

103

101

104

m=10000, ⌫=0.5

104

102

10 1 310

if k > kre + Ki and F (αk ) > F (αk−1 ) thenkre ← k, Ki+1 ← 2Ki , i ← i + 1# ‘mt’ηd ← δ · ηd + (1 − δ) · 1# ‘st’
k+1
k−1
k
k−1
β
←α
, α ←α# ‘re’
tk+1 ← 1, tk ← 0# ‘re’
end ifmin φU (w)

105

100

# Step 2’# Step 3∥w∥2 ≤1n=1000, ⌫=0.5
SeDuMiLIBSVM
FAPG
LIBLINEAR

Time (sec.)ϛχϚοΫε֬཰ϚγϯU ′ , x(α)()
U ′ := (α+ , α− ) ∈ R2n | ∥αo ∥2 ≤ κ, o ∈ {+, −}
%
&
%
&
1/2
1/2
x(α) = x̄+ + Σ+ α+ − x̄− + Σ− α−Time (sec.)

Ϟσϧ

(3)

2

統計的モデリングと計算アルゴリズムの数理と展開

報告書

空間の低ランク性と平滑性を考慮した
フーリエ係数最適化による MR 超解像
本谷秀堅†

河村 直輝†

横田 達也†† 名古屋工業大学 〒 466-8555 愛知県名古屋市昭和区御器所町
あらまし

MR 画像は撮影時に空間分解能を高く設定すると，S/N 比が著しく低下する．本稿では同一試料を複数方

向から撮影した MR 画像を超解像化することで，S/N 比を保持したまま高分解能な MR 画像を構成する．提案法では，対象の輪郭情報を用いる超解像法に TV 正則化と低ランク化を導入する．
キーワード 画像処理，超解像，Total Variation, 低ランクテンソル補完，ADMM

Hidekata HONTANI† , Naoki KAWAMURA† , and Tatsuya YOKOTA†† Nagoya Institute of Technology

1. は じ め に
磁気共鳴画像法 (MRI) は生体内部の 3 次元空間情報を高コントラストで観察するために有用であり，重要なモダリティの一つとして広く利用されている．MRI では体内の水素原子濃度を磁気を用いて観測し，体内断面をスライス画像として撮影する．スライス厚を小さくすれば空間分解能を高くすることが

図 1 提案法の流れ．MRI により対象物体は異なる複数方向からそれできるが，その場合はスライス毎の水素原子含有量が減少し，

ぞれ非等方的に観測される. 各画像を統合しても，観測されない

S/N 比が大きく低下する．そのため MRI では，十分な S/N 比

周波数領域があり，これら未知の高周波成分を超解像で復元する.

を確保するためにスライス厚を大きくとって撮影する．ただしスライス厚を大きくすると，MR 画像は磁場方向の空間分解能保ちつつ高解像な MR 画像を構築できることが期待できる．しが低く，非等方的な 3 次元画像となる．すなわち，MRI の磁場

かし，軸から離れた高周波成分は観測できない．そこで，未知方向の高周波成分が欠損した状態で画像が得られる．

の高周波成分を補完するために新たに超解像法を提案し，効果

空間分解能が非等方的な MR 画像を画像処理する場合，既存的な高解像度化を目指す．

のアルゴリズムをそのまま適用できないケースが多々ある．そ2. 問 題 設 定のため，MR 画像を等方的に高解像度化することは非常に重要である．
超解像は画像信号における未知の高周波成分を復元する技術である．Gerchberg は対象の輪郭情報が既知の場合に，効果的に解くことができる超解像法を提案した [2]．Gerchberg の手法では信号領域と周波数領域の双方から誤差エネルギーを反復的に除去することで高解像度化を行う．一方で近年，Total

Variation(TV) 正則化による超解像法が注目を浴びている．TV

本節では，提案手法について説明する．同一対象を異なる 3 つの方向から撮影した MR 画像を I1 ∈ RN ×N ×n , I2 ∈ Rn×N ×N ,I3 ∈ RN ×n×N とする．これらの磁場方向は互いに直交し，
n < N である．すなわち，各スライス画像サイズは N × N であり，スライス枚数は n である．I1 ∼ I3 から，等方的で高分解能な画像 I ∈ RN ×N ×N を復元したい．ここで，Ii のフーリエ変換を Fi とすると，目標の F 周波数空間のうち Fi の成分

正則化は画像のエッジを保持しつつ滑らかさを評価できるため，

は図 1 に示すように，磁場方向に向かって狭帯域に計測されて超解像に有効である [3]．また，最近では行列の低ランク化によいる．この観測領域を Ωi とする．3 つの観測領域を統合した

るテンソル補完技術も，その有用性が注目されている [4]．

Ω := Ω1 ∪ Ω2 ∪ Ω3 は F 周波数空間全てを表現できず，未知の提案法では，図 1 のように異なる方向から撮影した複数枚のMR 画像を用いる．互いに異なる磁場方向から撮影した数枚の高周波成分が存在する．この未知の高周波成分を，次節で提案する超解像により復元する．

MR 画像は，低分解能な方向を相互に補い合うため，S/N 比を—1—

berg 法に TV 正則化とテンソルの低ランク補完を組み合わせる．提案法では，次のような凸最適化問題を解くことで，画像
3

x = vec(I) ∈ RN ×1 を復元する．
f

j=1

+

(xj − Dj Hj x) + λT V ||x||T V + λLR

3
!
ρ
j=1

30

25

25

20

20

15
LRTVG
LRTV
TVG
TV
Gerchberg
ZP
sinc
NN

10

5

LRTVG
LRTV
TVG
TV
Gerchberg
ZP

5

sinc
NN

0
4

5

6

7

8

9

10

11

0

12

0.01

0.02

(a)

3

0.03

0.04

0.05

0.06

0.07

0.08

0.09

Noise Level [%]β3!
||M(j) ||tr
j=1

15

10

0

(b)
30

25

µ
||x − vec(Mj )||22 + ||f0 − RΩ f ||22 ,2
2

20

PSNR [dB]

min

3
!

MSE

我々の提案する超解像法では，対象の輪郭を用いる Gerch-

30

MSE

3. 提 案 法

15

10s.t. x = Gf , 0 = RΓ̄ x, M(j) = unfoldj (Mj ).

(1)

5

0
LRTVG LRTV

TVG

ここで，vec(·) はテンソルからベクトルへの変換，unfoldj (·) は

3 次元テンソルを方向 j に沿って行列に展開する．Mj はテン
ソル I の方向 j に対応する低ランク化のスラック変数であり，
2

f0 は f の観測された初期値である．xj = vec(Ij ) ∈ RN n×1

であり，Hj , Dj はそれぞれ磁場方向に向かって平滑化，サンTV Gerchberg ZP

sinc

NN

(c)
図 2

WEIZMANN デ ー タ セ ッ ト に よ る シ ミ ュ レ ー シ ョ ン 結 果 ．
(a)sharp image における MSE w.r.t. β ，(b)sharp image のノイズ強度変化における MSE(β = 4)，(c)β = 4，データセッ
ト集合で計測した PSNR.

プルを行う行列である．|| · ||TV は Total Variation，|| · ||tr は行列のトレースノルムであり，G は逆フーリエ基底行列であ33る．また RΓ̄ ∈ {0, 1}N ×N は画像中の対象領域外部を指定し，
3

3

RΩ ∈ [0, 1]N ×N は周波数空間の領域 Ω を指定する行列である．λT V , λLR , ρ, µ > 0 はそれぞれ TV 正則化，低ランク化，
スラック変数のフィッティング，周波数フィッティングの 4 つの項のバランスを調整するパラメータである．式 (1) は PDSや ADMM，MM アルゴリズム等で解く．ADMM では，式 (1)(a) 観測画像

を x, Mj , f と 2 つのラグランジュ係数について最適化する．

4. 実

験

WEIZMANN データセット [1] に対してシミュレーション実験を行った．N × N ，N = 120 の元画像を目標画像とし，これを 1 方向に向かって間隔 β でダウンサンプルした画像 I1 ，I2 ，(n = N/β) を得る．狭帯域な観測画像 I1 ，I2 から元画像を復元し，対象領域内部の誤差で従来法と比較評価を行った結果を
図 2 に示す．比較する手法は最近傍法 (NN)，Gerchberg 法 [2]，TV 正則化超解像法 (TV) [3]，低ランク補間を導入した TV 超
解像 (LRTV) [4]，提案法 (LRTVG)，提案法から低ランク補間項を除いたもの (TVG) である．提案法は従来法と比較して，βを変化させても精度の高い結果が得られた．しかし，強いノイズを付加すると従来法より性能が減少した．次にスダチの MR
実画像を用いて実験を行った結果を図 3 に示す．観測時には欠損していたエッジ成分が復元され，高解像度化された．

5. お わ り に
本稿では，同一試料を複数方向から撮影した MR 画像を超解像化することで，等方的に高分解能な MR 画像を構成した．

Gerchberg の超解像法を TV 正則化と低ランク化することで，効果的に超解像化できることを確認した．
謝辞

本研究は JSPS 科研費 26108003，及び 15K16067 の(b) 復元結果
図3

T2-MRI で撮影したスダチの実験結果.

文

献[1] S. Alpert, M. Galun, R. Basri, and A. Brandt, Image segmentation by probabilistic bottom-up aggregation and cue
integration, Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, vol. 0, pp. 1-8, 2007.
[2] R. W. Gerchberg. “Super-resolution through error energy
reduction.” Journal of Modern Optics, vol.21, no.9, pp709720, 1974.
[3] A. Marquina, and S. J. Osher. “Image super-resolution by
TV-regularization and Bregman iteration.” Journal of Scientific Computing, vol.37, no.3, pp.367-382, 2008.
[4] F. Shi, J. Cheng, L. Wang, P-T. Yap, and D. Shen. “LRTV:
MR image super-resolution with low-rank and total variation regularizations.” IEEE transactions on Medical Imaging, vol.34, no.12, pp.2459-2466, 2015.

助成を受けたものである．

—2—

େࣗ༝౓‫ܥ‬ͷσʔλಉԽͷͨΊͷ 2nd-order adjoint ๏Λ༻͍ͨߴ଎ෆ࣮֬ੑධՁ๏
ҏ౻৳Ұ 1 ɺ௕ඌେಓ 1,2 ɺՃೲকߦ 1 ɺࠇՏఱ 212౦‫ژ‬େֶ஍਒‫ॴڀݚ‬౦‫ژ‬େֶେֶӃ৘ใཧ޻ֶ‫ڀݚܥ‬ՊσʔλಉԽ [1] ͸γϛϡϨʔγϣϯϞσϧͱ‫؍‬ଌσʔλΛϕΠζ౷‫ֶܭ‬Λ༻͍ͯ༥߹͢Δ‫جࢉܭ‬൫ٕज़Ͱ͋Γɺࣄ‫ޙ‬෼෍ͷ‫ܗ‬ঢ়ΛධՁ͢Δ͜ͱʹΑΓɺ‫؍‬ଌෆՄೳͳ‫ܥ‬ͷ಺෦ঢ়ଶͷਪఆ΍ϞσϧύϥϝʔλͷਪఆΛߦͳ͏͚ͩͰͳ͘ɺͦͷਪఆ஋ͷෆ࣮֬ੑΛซͤͯධՁ͢Δ͜ͱ͕ՄೳͰ͋
Δɻࣄ‫ޙ‬෼෍ͷߏ੒ʹ༻͍ΒΕΔσʔλಉԽ͸ஞ࣍‫ͱܕ‬ඇஞ࣍‫ʹܕ‬େผ͞ΕɺΞϯαϯϒϧΧϧϚ
ϯϑΟϧλ [2] ΍ཻࢠϑΟϧλ [3] ౳ʹ‫ͮ͘ج‬ஞ࣍σʔλಉԽ͸ɺ͞·͟·ͳՊֶ෼໺Ͱ༻͍ΒΕ
͍ͯΔɻ͔͠͠ɺ͜ΕΒͷख๏͸ϞϯςΧϧϩϕʔεͷख๏Ͱ͋ΓҰൠʹਪఆ͓Αͼෆ࣮֬ੑධՁ
ʹඞཁͳ‫ྔࢉܭ‬ɾϝϞϦ͕γϛϡϨʔγϣϯϞσϧͷࣗ༝౓ͷࢦ਺ΦʔμʔͱͳΔͨΊɺγϛϡϨʔγϣϯϞσϧͷ‫ن‬໛͕େ͖͘ͳΔͱධՁ͕೉͘͠ͳΔɻҰํͰɺඇஞ࣍σʔλಉԽ๏Ͱ͋Δ̐࣍‫ݩ‬ม෼๏ [4] ͸ɺਪఆʹඞཁͳ‫ྔࢉܭ‬ɾϝϞϦΛϞσϧͷࣗ༝౓ͷઢ‫ܗ‬Φʔμʔʹ཈͑Δ͜ͱ͕Ͱ͖ΔͨΊɺେ‫ن‬໛γϛϡϨʔγϣϯϞσϧʹద͍ͯ͠Δɻ͔͠͠ɺैདྷͷ̐࣍‫ݩ‬ม෼๏ͷ࿮૊ΈͰ͸ޯ഑๏ʹΑΔࣄ‫࠷཰֬ޙ‬େղͷΈ͕෼͔Δ͚ͩͰ͋Γɺͦͷෆ࣮֬ੑΛධՁ͢Δ͜ͱ͸‫ݪ‬ཧతʹෆՄೳͰ͋ΔͨΊɺஞ࣍σʔλಉԽ๏ͱ૊Έ߹ΘͤΔ౳ͷ͞·͟·ͳΞυϗοΫͳ޻෉Λ‫ڽ‬Β͢͜ͱʹΑΓɺෆ࣮֬ੑͷධՁΛߦͳ͍ͬͯͨɻ
ͦ͜Ͱզʑ͸ɺͦͷ໰୊Λղܾ͢ΔͨΊʹɺ2nd-order adjoint ๏ [5] ͱ͍͏‫ࢉܭ‬๏Λར༻ͯ͠ɺ
େ‫ن‬໛ͳγϛϡϨʔγϣϯϞσϧʹରͯ͠΋ߴ଎͔ͭߴਫ਼౓ͳෆ࣮֬ੑධՁΛՄೳʹ͢ΔΞϧΰϦ
ζϜΛ։ൃͨ͠ [6]ɻෆ࣮֬ੑ͸ࣄ‫ޙ‬෼෍Λۙࣅ͢Δਖ਼‫ن‬෼෍ͷ‫ڞ‬෼ࢄߦྻͷର֯੒෼ɺͭ·Γϔοη‫ྻߦٯ‬ͷର֯੒෼Ͱۙࣅ͞ΕΔ͕ɺຊ‫Ͱڀݚ‬͸ɺ࣮༻্ෆ࣮֬ੑΛධՁ͢΂͖ม਺ͷ਺͕γϛϡϨʔγϣϯϞσϧͷࣗ༝౓ʹൺ΂ͯѹ౗తʹগͳ͍͜ͱʹ஫໨͠ɺ2nd-order adjoint ๏Λ͏·͘ར༻͢Δ͜ͱͰɺ஫໨͢Δཁૉͷෆ࣮֬ੑͷΈΛ௚઀‫ٻ‬ΊΔΞϧΰϦζϜΛߏஙͨ͠ɻ2nd-orderadjoint ๏Λ༻͍ͳ͍ैདྷ·Ͱͷෆ࣮֬ੑͷ‫ࢉܭ‬๏Ͱ͸ඞཁͱ͢ΔϝϞϦ͕γϛϡϨʔγϣϯϞσ
ϧͷࣗ༝౓ͷ̎৐ɺ‫͕ྔࢉܭ‬γϛϡϨʔγϣϯϞσϧͷࣗ༝౓ͷ̏৐Ҏ্ඞཁͱ͢Δͷʹରͯ͠ɺఏҊख๏͸ϝϞϦ͕γϛϡϨʔγϣϯϞσϧͷࣗ༝౓ͷ̍৐ɺ‫͕ྔࢉܭ‬γϛϡϨʔγϣϯϞσϧͷ
ࣗ༝౓ʹґΒͳ͍ (γϛϡϨʔγϣϯϞσϧͷ࣌‫ྻܥ‬Λ‫͢ࢉܭ‬Δ‫ͱྔࢉܭ‬ಉఔ౓) ʹͳΔͷͰɺେ‫ن‬໛ͳγϛϡϨʔγϣϯϞσϧʹରͯ͠ඇৗʹ༗རͰ͋ΔɻఏҊख๏ͷ‫ূݕ‬ͷͨΊɺ‫ݻ‬ମ-ӷମ૬ք໘ͷҠಈμΠφϛΫεΛ‫ه‬ड़͢ΔϑΣʔζϑΟʔϧυϞ

1σϧͷ࣌ؒൃల͔ΒಘΒΕΔٖࣅσʔλʹຊఏҊख๏Λద༻͠ɺϞσϧͷॳ‫ظ‬ঢ়ଶ͓ΑͼϞσϧύ
ϥϝʔλͷਪఆɺϞσϧύϥϝʔλͷෆ࣮֬ੑධՁΛߦͳͬͨɻ‫ূݕ‬ͷ݁ՌɺຊఏҊख๏͸Ծఆ͞Εͨॳ‫ظ‬ঢ়ଶ͓ΑͼϞσϧύϥϝʔλΛਖ਼͘͠ਪఆ͢Δ͜ͱ͕Ͱ͖ɺ͞Βʹɺσʔλͷྔɾ࣭ʹԠͨ͡ෆ࣮֬ੑΛఆྔత͔ͭߴ଎ʹධՁ͢Δ͜ͱʹ੒ޭͨ͠ɻຊఏҊख๏ʹΑͬͯಘΒΕΔෆ࣮֬ੑ͸ɺ࣮‫ݧ‬౳΁ͷϑΟʔυόοΫʹର͢Δॏཁͳ஌‫ݟ‬Λ༩͑Δ͜ͱ͕Ͱ͖ɺྫ͑͹ɺඞཁͳਫ਼౓ͷਪఆ݁ՌΛಘΔʹ͸Ͳͷఔ౓σʔλ͕ඞཁ͔ͳͲΛ‫ͯ͑͘ڭ‬ΕΔɻখ‫ن‬໛ࣗ༝౓ͷγϛϡϨʔγϣϯϞσϧʹରͯ͠͸ஞ࣍σʔλಉԽ๏౳Λ࢖͏͜ͱͰͦΕ͸Մೳ͕ͩͬͨɺେ‫ن‬໛ࣗ༝౓ͷγϛϡϨʔγϣϯϞσϧʹରͯ͠΋ɺຊఏҊख๏ʹΑΓͦΕ͕Մೳʹͳͬͨɻࢀߟจ‫ݙ‬[1] S. Reich and C. Cotter, Probabilistic Forecasting and Bayesian Data Assimilation (Cambridge University Press, Cambridge, 2015).
[2] G. Evensen, The Ensemble Kalman Filter: theoretical formulation and practical implementation, Ocean Dynamics 53, 343 (2003).
[3] G. Kitagawa, Introduction to Time Series Modeling, Chapman & Hall/CRC Monographson Statistics & Applied Probability (CRC Press, Boca Raton, FL, 2010).
[4] F.-X. Le Dimet and O. Talagrand, Variational algorithms for analysis and assimilation of
meteorological observations: theoretical aspects, Tellus A 38A, 97 (1986).
[5] F.-X. Le Dimet, I. M. Navon, and D. N. Daescu, Second-order information in data assimilation, Monthly Weather Review 130, 629 (2002).
[6] S. Ito, H. Nagao, A. Yamanaka, Y. Tsukada, T. Koyama, M. Kano, and J. Inoue, Dataassimilation for massive autonomous systems based on a second-order adjoint method,
Physical Review E 94, 043307 (2016).

2

座標ごとの変数変換によって得られる Stein 型の分布
清 智也 （東大 情報理工）∗概要与えられた多次元の確率分布に対し，座標ごとの確率積分変換によってコピュラが得られる
ことはよく知られている（Sklar の定理）。本研究では，座標ごとの変数変換によって，Stein型の等式を満たす分布が一意的に得られることを示す。ただしいくつかの正則条件を要する。
本結果は Marshall and Olkin (1968, Numer. Math.) による対角スケーリング定理の非線形版と解釈できる。

1 Stein 型の等式
d を正の整数とする。Rd 上の確率分布 µ のうち，各周辺分布 µi (i = 1, . . . , d) が絶対連続で，!
x dµi = 0 かつ R x2i dµi < ∞ となるようなものの全体を P 2 とおく。µi は絶対連続だが µ 自R i!

身は絶対連続とは限らないことに注意する。また，局所絶対連続な関数 f : R → R のうち，その導関数 f ′ が本質的に有界であるようなものの全体を B とおく。

定義 1. 分布 µ ∈ P 2 が Stein 型分布であるとは，全ての f ∈ B に対して

"

⎛
⎞
"
d
%
⎝
⎠
f (xi )
xj dµ = f ′ (xi )dµ,i = 1, . . . , d,

(1)

j=1

を満たすこととする。また式 (1) を Stein 型の等式と呼ぶ。

!
f (x1 )x1 dµ = f ′ (x1 )dµ に帰着され，µ は標(d
準正規分布に限られる（例えば [1]）。同様に，独立な分布 µ = i=1 µi で Stein 型となるのは各
d = 1 の場合，式 (1) はいわゆる Stein の等式!

µi が標準正規分布のときのみである。多変量正規分布 µ = N(0, S) の場合，µ が Stein 型となるの)d
は各 i に対して j=1 Sij = 1 のときであることが示される。
もし確率変数ベクトル (X1 , . . . , Xd ) が Stein 型分布に従うならば，式 (1) より，それらの和)j Xj は各 Xi の任意の単調増加変換 f (Xi ) と正の相関を持つことが分かる。この性質の応用として，多変量データを総合指数にまとめる方法が議論されている [4]。
座標ごとの変換T (x) = (T1 (x1 ), . . . , Td (xd )),x = (x1 , . . . , xd ) ∈ Rd ,∗ 本研究は科研費 (課題番号: 26108003, 26540013) の助成を受けたものである。

のうち，各 Ti : R → R ∪ {±∞} が非減少関数で，かつ像測度 T♯ µ が P 2 に属すようなものの全体
を Tcw (µ) とおく。ここで cw は coordinate-wise の略である。

座標ごとの変換で得られる分布としてコピュラがよく知られている [3]。本研究では，与えられた µ ∈ P 2 に対し，T♯ µ が Stein 型となるような座標ごとの変換 T ∈ Tcw (µ) を見つける問題を考える。そのような T を µ の Stein 型変換と呼ぶことにしよう。

2 主結果
定理 1. µ ∈ P 2 とし，T ∈ Tcw (µ) とする。このとき次の 2 つの条件は同値である：(i) T は µ の Stein 型変換である。
(ii) T♯ µ は集合 {S♯ µ | S ∈ Tcw (µ)} において E を最小化する。
定理 2 (一意性). µ ∈ P 2 のサポートが周辺分布 µ1 , . . . , µd のサポートの直積に一致すると仮定する。このとき µ の Stein 型変換は高々一つである。
分布 µ ∈ P 2 が共正値 (copositive) であるとは，が成り立つこととする。

! )(T )2 dµ) !i i 2
inf> 0.
T ∈Tcw (µ)
Ti dµi

定理 3 (存在性). µ が共正値ならば，µ の Stein 型変換が存在する。
定理 2，定理 3 は，次の対角スケーリング定理の一般化とみなすことができる。
定理 4 (Marshall and Olkin [2]). S は d 次の半正定値行列で，かつ真に共正値 (strictly copositive)inf

w1 ,...,wd >0

) )
i

)j

wi Sij wj

2
i wi

>0

であるとする。このとき，ある正の数 w1 , . . . , wd で，各 i についてのが一意的に存在する。

)d

j=1 wi Sij wj = 1 となるも

参考文献
[1] Chen, L. H. Y., Goldstein, L., and Shao, Q. (2011). Normal Approximation by Stein’s Method,
Springer.
[2] Marshall, A. W., Olkin, I., (1968). Scaling of matrices to achieve specified row and column sums.
Numer. Math., 12, 83–90.
[3] Nelsen, R. B. (2006). An Introduction to Copulas, 2nd ed., Springer.
[4] Sei, T. (2016). An objective general index for multivariate ordered data, J. Multivariate Anal.,
147, 247–264.

An operational characterization of the notion of probability byalgorithmic randomness and its applications⇤
Kohtaro Tadaki
Department of Computer Science, College of Engineering, Chubu University
1200 Matsumoto-cho, Kasugai-shi, Aichi 487-8501, Japan
E-mail: tadaki@cs.chubu.ac.jp
http://www2.odn.ne.jp/tadaki/

The notion of probability plays an important role in almost all areas of science and technology. In modern mathematics, however, probability theory means nothing other than measure
theory [10], and the operational characterization of the notion of probability is not establishedyet. In this talk, based on the toolkit of algorithmic randomness we present an operational
characterization of the notion of probability, called an ensemble.
Algorithmic randomness, also known as algorithmic information theory, is a field of mathematics which enables us to consider the randomness of an individual infinite sequence [14, 11, 2,
12, 3, 4, 13, 7]. We use the notion of Martin-Löf randomness with respect to Bernoulli measure[12] to present the operational characterization. As the first step of the research of this line, inthis talk we consider the case of finite probability space, i.e., the case where the sample spaceof the underlying probability space is finite, for simplicity.
We give a natural operational characterization of the notion of conditional probability interms of ensemble, and give equivalent characterizations of the notion of independence between
two events based on it. Furthermore, we give equivalent characterizations of the notion of
independence of an arbitrary number of events/random variables in terms of ensembles. In
particular, we show that the independence between events/random variables is equivalent to the
independence in the sense of van Lambalgen’s Theorem [17], in the case where the underlying
finite probability space is computable.
In the talk we make applications of our framework to quantum mechanics, information
theory, and cryptography in order to demonstrate the wide applicability of our framework to
the general areas of science and technology.
For the detail of this talk, see Tadaki [15].
Key words: probability, algorithmic randomness, operational characterization, Martin-Löfrandomness, Bernoulli measure, conditional probability, independence, van Lambalgen’s Theorem, quantum mechanics, information theory, cryptography

References[1] R. B. Ash, Information Theory. Dover Publications, Inc., New York, 1990.
[2] G. J. Chaitin, “On the length of programs for computing finite binary sequences,” J. Assoc.
Comput. Mach., vol. 13, pp. 547–569, 1966.
⇤
This work was partially supported by JSPS KAKENHI Grant Numbers 24540142, 15K04981. This work
was partially done while the author was visiting the Institute for Mathematical Sciences, National University ofSingapore in 2014.

1[3] G. J. Chaitin, “A theory of program size formally identical to information theory,” J.
Assoc. Comput. Mach., vol. 22, pp. 329–340, 1975.
[4] G. J. Chaitin, Algorithmic Information Theory. Cambridge University Press, Cambridge,
1987.
[5] B. S. DeWitt and N. Graham (eds.), The Many-Worlds Interpretation of Quantum Mechanics. Princeton University Press, Prenceton, 1973.
[6] P. A. M. Dirac, The Principles of Quantum Mechanics, 4th ed. Oxford University Press,
London, 1958.
[7] R. G. Downey and D. R. Hirschfeldt, Algorithmic Randomness and Complexity. SpringerVerlag, New York, 2010.
[8] H. Everett, III, ““Relative State” formulation of quantum mechanics,” Rev. Mod. Phys.,vol. 29, no. 3, pp. 454–462, 1957.
[9] J. Katz and Y. Lindell, Introduction to Modern Cryptography. Chapman & Hall/CRC Press,
2007.
[10] A. N. Kolmogorov, Foundations of the theory of probability, Chelsea Publishing Company,
New York, 1950.
[11] A. N. Kolmogorov, “Three approaches to the quantitative definition of information,” Problems Inform. Transmission, vol. 1, no. 1, pp. 1–7, 1965.
[12] P. Martin-Löf, “The definition of random sequences,” Inform. and Control, vol. 9, pp. 602–
619, 1966.
[13] A. Nies, Computability and Randomness. Oxford University Press, Inc., New York, 2009.
[14] R. J. Solomono↵, “A formal theory of inductive inference. Part I and Part II,” Inform. andControl, vol. 7, pp. 1–22, 1964; vol. 7, pp. 224–254, 1964.
[15] K. Tadaki, “An operational characterization of the notion of probability by algorithmic
randomness and its applications,” arXiv:1611.06201, November 2016. Available at URL:
https://arxiv.org/abs/1611.06201[16] K. Tadaki, “A refinement of quantum mechanics by algorithmic randomness and its application to quantum cryptography,” Proceedings of the 2017 Symposium on Cryptography
and Information Security (SCIS2017), 4A1-2, January 24-27, 2017, Naha, Japan.
[17] M. van Lambalgen, Random Sequences. Ph.D. dissertation, University of Amsterdam, 1987.
[18] R. von Mises, Probability, Statistics and Truth. Dover Publications, Inc., New York, 1957.
[19] R. von Mises, Mathematical Theory of Probability and Statistics. Academic Press Inc., New
York, 1964.

2ʦใࠂॻʧ
௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍ͷࣹӨΞϧΰϦζϜʹ
ΑΔ୳ࡧ͍ͭͯ
On the Search by the Projection Algorithm for the Output Distributionthat Achieves the Channel Capacityத઒ ݈࣏ɼ౉෦ ߁ฏɼࠤ෢ ୓ే௕Ԭٕज़Պֶେֶ

1. ֓ཁຊ࿦จͰ͸཭ࢄతແ‫ه‬Ա௨৴࿏ͷ௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍Λ୳ࡧ͢ΔΞϧΰϦζϜʹ
͍ͭͯߟ࡯ͨ͠ɻ͜ͷΞϧΰϦζϜ͸ɼग़ྗ෼෍શମͷۭؒʹ͓͚Δ̍ͭͷग़ྗ෼෍͔ΒΞϑΟϯ෦෼ۭؒ΁ͷࣹӨͷ൓෮Ͱ͋Δɻ௨৴࿏༰ྔͷ໰୊͸ɼϢʔΫϦουۭؒʹ͓͚Δ༗‫ݸݶ‬ͷ఺ʹର͢Δ࠷খแ‫ؚ‬ԁ໰୊ͱಉ༷ͷ‫ز‬Կֶతߏ଄Λ༗͢ΔɻϢʔΫϦουۭؒ಺ͷई౓͸ϢʔΫϦου‫͋Ͱ཭ڑ‬Γɼग़ྗ෼෍શମͷۭؒͰͷई౓͸ Kullback-Leibler ৘ใྔͰ͋Δɻຊ࿦จͰ͸ɼ·ͣɼϢʔΫϦουۭؒʹ͓͚Δ࠷খแ‫ؚ‬ԁ໰୊Λߟ࡯͠ɼ࠷খแ‫ؚ‬ԁͷத৺Λ‫͢ࢉܭ‬ΔΞϧΰϦζϜΛ։ൃͨ͠ɻͦͯ͠ɼͦ͜ͰಘΒΕͨΞϧΰϦζϜΛ௨৴࿏༰ྔͷ໰୊ʹద༻͠ɼ௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍Λ‫͢ࢉܭ‬ΔΞϧΰϦζϜΛಘͨɻ཭ࢄతແ‫ه‬Ա௨৴࿏ͷ௨৴࿏༰ྔ C ͸૬‫ޓ‬৘ใྔͷ࠷େ஋ͱͯ͠ఆٛ͞ΕΔɻC ͸·ͨɼ

Kullback-Leibler (KL) ৘ใྔʹؔ͢Δ min max ໰୊ͷղͱͯ͠΋ఆࣜԽ͞ΕΔɻ௨৴࿏༰ྔʹؔ͢Δ min max ໰୊ʹ͓͚Δ KL ৘ใྔΛϢʔΫϦου‫ʹ཭ڑ‬ஔ͖‫͑׵‬ΔͱɼϢʔΫϦουۭؒʹ͓͚Δྨࣅͷ໰୊͕ಘΒΕΔɻͦΕ͸ɼ༗‫ʹ߹ू఺ݶ‬ର͢Δ࠷খแ‫ؚ‬ԁͷ໰୊Ͱ͋Δɻ
௨৴࿏༰ྔͷ໰୊͸ग़ྗ෼෍ʹؔ͢Δ࠷దԽ໰୊Ͱ͋Γɼߟ͑Δۭؒ͸ϢʔΫϦουۭؒͰ͸ͳ֬͘཰෼෍શମ͕ͳ͢ू߹Ͱ͋Δɻͦͷू߹্ͷ‫ز‬Կֶ͸৘ใ‫ز‬ԿͰ͋ΔɻϢʔΫϦου‫ز‬Կͱ৘ใ‫ز‬Կ͸‫؁‬རͷ α ‫ز‬Կͷཱ৔͔Β౷Ұతʹߟ͑Δ͜ͱ͕Ͱ͖Δɻ͕ͨͬͯ͠ɼ࠷খแ
‫ؚ‬ԁ໰୊ͷղ๏ΞϧΰϦζϜΛɼ‫ز‬ԿతͳྨࣅੑΛ௨ͯ͠ɼ௨৴࿏༰ྔ໰୊ʹద༻Ͱ͖Δɻͦͷࡍʹॏཁͳ͜ͱ͸ɼ࠷খแ‫ؚ‬ԁ໰୊Λߟ࡯͢ΔͷʹϢʔΫϦου‫ز‬ԿͷͲΜͳੑ࣭Λ࢖ͬͯ΋Α͍ͱ͍͏͜ͱͰ͸ͳ͘ɼ௨৴࿏༰ྔ໰୊ʹద༻Մೳͳੑ࣭͚ͩΛ࢖͏ඞཁ͕͋Δɼͱ͍͏
͜ͱͰ͋Δɻͦͷ݁Ռɼ࠷খแ‫ؚ‬ԁ໰୊ͰಘΒΕͨΞϧΰϦζϜΛ΄΅ࣗಈతʹ௨৴࿏༰ྔ໰
୊ʹҠ২͢Δ͜ͱ͕Ͱ͖Δɻຊ࿦จͰ͸ɼ࣮ࡍʹ྆ํͷ‫ز‬Կʹ‫ڞ‬௨ͷੑ࣭ͱͯ͠ɼ‫཭ڑ‬ɼॏ৺࠲ඪɼ಺ੵɼϐλΰϥεͷఆཧɼΞϑΝΠϯ෦෼ۭؒ΁ͷࣹӨͷΈΛ࢖ͬͯΞϧΰϦζϜΛ։ൃͨ͠ɻ͜ͷΞϧΰϦζϜΛຊ࿦จͰ͸ʮࣹӨΞϧΰϦζϜʯͱ‫ݺ‬Μͩɻ͞Βʹɼैདྷѻ͍͕ࠔ೉ͱ͞Ε͍ͯͨೖྗΞϧϑΝϕοτ਺͕ग़ྗΞϧϑΝϕοτ਺ΑΓଟ͍৔߹ʹ͓͍ͯɼͦͷࠔ೉͞Λ؇࿨ͯ͠ɼ௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍Λ‫͢ࢉܭ‬Δํ๏ΛఏҊͨ͠ɻ௨৴࿏༰ྔͷ໰୊ͷ‫ج‬൫ͱͳΔ৘ใ‫ز‬Կ͕ϢʔΫϦου‫ز‬ԿͱҟͳΔ఺͸ɼϢʔΫϦου‫ز‬ԿͰ͸ 1 ͭͷ࠲ඪ‫ܥ‬Λ༻͍Δͷʹରͯ͠ɼ৘ใ‫ز‬ԿͰ͸‫ʹ͍ޓ‬૒ରͳ 2 ͭͷ࠲ඪ‫ܥ‬Λ༻͍Δ͜
ͱͰ͋ΔɻϢʔΫϦου‫ز‬Կͱ৘ใ‫ز‬Կ͸ɼ࣮਺ α ʹΑͬͯ α ‫ز‬Կͱ͍͏౷߹తͳ‫ز‬Կߏ଄ͷ଒ͷҰ෦ͱͯ͠ߟ͑ΒΕΔɻϢʔΫϦου‫ز‬Կ͸ α = 0, ৘ใ‫ز‬Կ͸ α = ±1 ͱ͍͏ಛผͳ஋ʹରԠ͢Δ‫ز‬Կͱͯ͠ଊ͑Δ͜ͱ͕Ͱ͖Δɻα ‫ز‬ԿͰ͸ɼα divergence, ಺ੵɼPythagoras ͷఆཧɼα ࣹӨ͕࢖͑Δɻ࠷খแ‫ؚ‬ԁ໰୊ʹؔ͢Δఆཧͷূ໌ʹ͓͍ͯར༻ͨ͠ϢʔΫϦου‫ز‬Կͷੑ࣭͕‫཭ڑ‬ɼ಺ੵɼPythagoras ͷఆཧɼࣹӨ͚ͩͳͷͰɼͦ͜ͰಘΒΕͨఆཧ΍ΞϧΰϦζ
Ϝ͸༰қʹ௨৴࿏༰ྔͷ໰୊ʹద༻͢Δ͜ͱ͕Ͱ͖ͨɻ
ϢʔΫϦου‫ز‬Կ͸զʑʹͳ͡Έͷ‫ز‬ԿͳͷͰɼΞϧΰϦζϜͷ։ൃͱͦͷূ໌͕͠΍͍͢ɻͦΕΒ͸͋Δఔ౓‫ڧ‬Ҿͳํ๏Ͱ΋Ͱ͖ͯ͠·͏͕ɼͦΕ͸͋ΔҙຯͰ͍͍͜ͱͰ͋Δɻ͔͠͠ɼ
ͦΕͱಉ͡ํ๏Λ৘ใ‫ز‬ԿʹҠ২͠Α͏ͱ͢Δͱɼ͏·͍͔͘ͳ͍৔߹͕͋ͬͨɻͦ͜Ͱɼͦ
Ε͸‫ڧ‬ҾͰࣗવͳํ๏Ͱ͸ͳ͔ͬͨͱ‫͕ͭ͘ؾ‬ɻ৘ใ‫ز‬Կʹద༻Ͱ͖ΔΑ͏ʹߟ͑௚͢ͱɼ‫ٯ‬ʹϢʔΫϦου‫ز‬Կʹؐ‫͖ͯͰݩ‬ɼϢʔΫϦου‫ز‬ԿͰผͷํ๏Λൃ‫͢ݟ‬Δ͜ͱ͕Ͱ͖ͨɻͭ
·ΓɼϢʔΫϦου‫ز‬Կͱ৘ใ‫ز‬Կʹ‫ڞ‬௨ͷੑ࣭͚ͩΛ࢖ͬͨΞϧΰϦζϜ͕࠷΋ࣗવͰ͋Γɼ࣮͸࠷΋؆୯ͳํ๏Ͱ͋Δɻ͜Ε͸ඇৗʹ‫ڵ‬ຯਂ͍݁ՌͰ͋ΔɻͨͿΜɼॳΊ͔Β௨৴࿏༰ྔ
ͷ໰୊Λ৘ใ‫ز‬Կ͚ͩͰߟ͍͑ͯͨͷͰ͸ຊ࿦จͷࣹӨΞϧΰϦζϜ͸ಘΒΕͳ͔ͬͨͱࢥ͏ɻ
ϢʔΫϦου‫ز‬Կͱ৘ใ‫ز‬Կ͕༗‫ػ‬తʹ݁ͼ͍ͭͨ݁Ռͩͱࢥ͏ɻ

2. ຊ࿦จͰಘΒΕͨ݁Ռຊ࿦จͰಘΒΕͨ݁Ռ͸Լ‫ه‬ͷ௨ΓͰ͋Δɻ

• ௨৴࿏༰ྔͷ໰୊͸ϢʔΫϦουۭؒ Rn ʹ͓͚Δ࠷খแ‫ؚ‬ԁ໰୊ͱྨࣅͨ͠‫ز‬Կߏ଄Λ࣋ͭ͜ͱΛ໌Β͔ʹͨ͠ɻͦͷྨࣅੑʹ‫͍ͯͮج‬௨৴࿏༰ྔΛ‫͢ࢉܭ‬ΔͨΊͷࣹӨΞϧΰϦζϜΛ։ൃͨ͠ɻ

• ௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍ΛࣹӨΞϧΰϦζϜʹΑͬͯ୳ࡧ͢ΔͨΊͷఆཧΛࣔ͠ɼͦͷఆཧʹ͓͍ͯॏ৺࠲ඪ͕ॏཁͳ໾ׂΛ࣋ͭ͜ͱΛ໌Β͔ʹͨ͠ɻ• ௨৴࿏ߦྻ Φ ͷߦϕΫτϧ͕ඞͣ͠΋ҰൠͷҐஔʹ͋Δͱ͸‫ݶ‬Βͳ͍৔߹ʹ΋ɼ৽ͨʹ! Λఆٛͯ͠ɼͦͷߦϕΫτϧ͕ҰൠͷҐஔʹ͋ΔΑ͏ʹͯ͠ɼࣹΦ ʹ͍ۙ௨৴࿏ߦྻ Φ! ͷ֤ߦϕΫτϧ͸ɼ΋ͱͷ Φ ͷ֤ߦϕΫτϧͷ࣍‫ݩ‬Λ্ӨΞϧΰϦζϜΛద༻ͨ͠ɻΦ͛ͨ΋ͷͰ͋Δɻ͠͹͠͹ɼ͋Δ໰୊ʹ͓͍ͯ࣍‫ݩ‬Λ্͛Δ͜ͱʹΑͬͯ໰୊ͷࠔ೉͕͞؇࿨͞ΕΔ৔߹͕͋Δɻຊ࿦จͰఏҊ͍ͯ͠Δํ๏͸࣍‫ݩ‬Λ্͛Δํ๏ͷͻͱͭͷྫͰ
͋Δɻͳ͓ɼຊ࿦จͷ੒ՌΛ·ͱΊͯ࿦จͱͯ͠ IEEE Transactions on Information Theory ʹ౤ߘͯ͠࠾࿥ͱͳΓɼ‫͞ࡌܝ‬ΕΔ༧ఆͰ͋Δɻ

3. ࠓ‫ޙ‬ͷ՝୊ࠓ‫ޙ‬ͷ՝୊ͱͯ͠Լ‫ه‬ͷ͜ͱ͕‫͛ڍ‬ΒΕΔɻ• ೚ҙͷೖྗΞϧϑΝϕοτ਺ m ʹରͯ͠௨৴࿏༰ྔΛୡ੒͢Δग़ྗ෼෍Λ‫͢ࢉܭ‬ΔΞϧΰϦζϜΛ࡞Δɻ• Ϩʔτ࿪Έؔ਺΍ capacity constraint function Λ‫͢ࢉܭ‬ΔࣹӨΞϧΰϦζϜΛ࡞Δɻ
• ௨৴࿏༰ྔΛ‫͢ࢉܭ‬Δஞ࣍ΞϧΰϦζϜͰ͋Δ༗ຊΞϧΰϦζϜΛ࠷খแ‫ؚ‬ԁͷ໰୊ʹҠ২͢Δɻ␨ ⍬␄ ໵⏀๲໱⎩ λ౫шߏ‫ڤ‬ഇ೘୔രλౠ ∗The rate-distortion (RD) theory studies and characterizes the performance oflossy compression systems by the RD function [1]. There are studies interpreting
learning and prediction problems by rate-distortion theoretic views. In particular,
clustering methods can be considered as vector quantizers and naturally related to
the RD theory. RD theoretic interpretations have been obtained also for problems
such as classification and sequence prediction.
In this study, we provide a RD theoretic view of Bayesian learning. We formulate
a rate-distortion problem by the distortion measure defined by the pointwise regretof the model parameter θ,p(xn |θ̂(xn )),d(x , θ) = logp(xn |θ)nwhere xn is the data set and θ̂ is the maximum likelihood estimator. We show thatthe generalized Bayesian posterior distribution appears as the optimal solution to
the problem. The generalized posterior distribution has been used for purposessuch as the posterior consistency [2] and model selection [3].
We also discuss the connection of this view to the asymptotic theory of Bayesianlearning, which characterizes the generalization error of a learning machine bya constant called learning coeﬃcient [4]. It is demonstrated that the learning
coeﬃcient provides an upper bound of the RD dimension [5], which is defined by
the asymptotic behavior of the rate-distortion function. Let λ be the learning
coeﬃcient. It is known that λ is the leading term of the negative log-marginal
likelihood, and if the model is regular, 2λ = m, the dimension of the parameter,whereas if it is non-regular, 2λ ≤ m holds in general [4]. The asymptotic behaviorof the RD function of the above problem shows that twice the learning coeﬃcient
provides an upper bound of the RD dimension,
dimRD ≤ 2λ.
which is defined by the RD function R(D) as dimRD = limD→0 −R(D)
1
log D
2

∗

൦ֹՠࡦѕӁ੒Ӂ: wkazuho@cs.tut.ac.jp

1

We also provide a RD theoretic interpretation of the Dirichlet process (DP)means clustering, which is a simple extension of the K-means clustering and estimates the number of clusters from data [6]. More specifically, DP-means runs theusual K-means with K = 1, and generates a new cluster when the distance from adata point to its nearest cluster center is larger than a constant η, which is called
the penalty parameter and prespecified by the user.
From an RD theoretic point of view, the penalty parameter can be considered
as controlling the maximum distortion in the training data set {x1 , · · · , xN }(xi ∈RL ) [7],max d(xi , θc(i) ),iwhere c(i) denotes the cluster label of the ith data point and θ1 , · · · , θK̂ are thecluster centers obtained by the algorithm. The logarithm of the estimated numberof clusters per dimension of data, log K̂/L, can be considered as the rate. The RD
theory of the maximum distortion criterion [8] implies thatlog K̂
→ R(η)
L
as N → ∞ and L → ∞, where R is the RD function of the distortion measureused in the algorithm [7].

References[1] T. Berger, Rate Distortion Theory: A Mathematical Basis for Data Compression, Englewood Cliﬀs, NJ: Prentice-Hall, 1971.
[2] T. Zhang, “Information theoretical upper and lower bounds for statisticalestimation,” IEEE Trans. on Inform. Theory, vol. 52, no. 4, pp. 1307–1321,
2006.
[3] M. Aitkin, “Posterior Bayes factors,” Journal of the Royal Statistical Society.
Series B, vol. 53, no. 1, pp. 111–142, 1991.
[4] S. Watanabe, Algebraic Geometry and Statistical Learning Theory, Cambridge
University Press, 2009.
[5] T. Kawabata and A. Dembo, “The rate-distortion dimension of sets and measures,” IEEE Trans. on Inform. Theory, vol. 40, no. 5, pp. 1564–1572, 1994.
[6] B. Kulis and M. I. Jordan, “Revisiting k-means: new algorithms via Bayesiannonparametrics,” in Proc. of ICML, 2012, pp. 513–520.
[7] M. Kobayashi and K. Watanabe, “A rate-distortion theoretic view of Dirichlet
process means clustering,” IEICE Trans. on Fundamentals, submitted.
[8] T. S. Han, Information-Spectrum Methods in Information Theory, Springer
Berlin Heidelberg, 2003.

2