Fused lasso に基づく非線形回帰モデリング1

九州大学大学院数理学府立石 正平

中央大学理工学部

小西 貞則

はじめに回帰分析において，複雑な構造を有する現象を解明するために，スプライン関数や動径基底関数に基づく非線形回帰モデルが用いられる．しかし，データを発生させた真の構造が滑らかでない場合，非線形回帰モデルの当てはめが困難となり，データから有効な情報を得ることが難しくなる．
本報告では，fused lasso (Tibshirani et al., 2005) を用いて構成した基底関数に基づく非線形回帰モデルを適用することで，突発的な変化をもつ非線形構造を推定する方法を提案する．基底関数展開に基づく非線形回帰モデルを仮定する場合，基底関数の調整パラメータの選択が本質的となる．この問題に対して，モデル評価基準を導出し，調整パラメータの選択を行い，提案手法の有効性を数値実験を通して検証する．

2

基底関数展開に基づく非線形回帰モデル
目的変数 Y と説明変数 X に関して，n 組のデータ {(yα , xα ); α = 1, 2, · · · , n} が観測されたとする．このとき，説明変数と目的変数の間に次の基底関数展開に基づくモデルを仮定する．

yα = wT b(xα ) + α ,(α = 1, 2, · · · , n) .

ここで，b(x)=(1,b1 (x),b2 (x),· · · ,bm (x))T は (m + 1) 次元基底関数ベクトル，w = (w0 , w1 ,w2 , · · · , wm )T は (m + 1) 次元係数パラメータベクトルとする．ここで，基底関数として次の式で与えられるガウス型基底関数を用いる.
{}(x − cj )2
bj (x) = exp −, j = 1, 2, · · · , m.
2νh2j(1)基底関数の広がりを調整パラメータ ν(> 0) によって調整することで滑らかな推定曲線が得られることから，その有効性が報告されている (Ando and Konishi, 2008)．

1ガウス型基底関数に含まれる中心ベクトル µ と関数の広がりの程度を表す h2 を決定する必要があり，k−means 法が用いられる場合がある．すなわち，n 個のデータ {x1 , · · · , xn }を，基底関数の個数である m 個のクラスタ C1 , C2 , · · · , Cm に分割し，各クラスタ Cj に含まれる nj 個のデータについて µj と h2j を次のように決定する．

ĉj =1 ∑xα ,nj x ∈Cαĥ2j =1 ∑(xα − cj )2 .
nj x ∈C
αj(2)j

しかし，このように構成された基底関数に基づく非線形回帰モデルは，滑らかな曲線になってしまい，平均構造が突然変化するデータに対する当てはめが困難となる．この問題を克服するために，本報告では，fused lasso の特殊な形である fused lasso signal approximator(FLSA; Friedman et al., 2007) を用いて基底関数を構成することを提案する．

3FLSA とガウス型基底関数の構築
n 組のデータ {(yα , xα ); α = 1, 2, · · · , n} に対して，FLSA は次の目的関数の最小化を行う．
nnn−1∑∑1∑2(yα − βα ) + λ1|βα | + λ2|βα − βα+1 |, α = 1, · · · , n.
2 α=1α=1α=1(3)

ここで， βα は，目的変数 yα の推定量であり，λ1 は β を真に 0 に推定するための正則化パラメータである．λ2 は隣り合う β の違いを真に 0 に推定するための正則化パラメータであり，λ2 の値を十分大きくすることによって, xFj = {xα |α ∈ Fj } に関する β は全て β̂FjFと推定される (∪nj=1
Fj = {1, 2, · · · , n}, Fk ∩ Fl = ∅ (k 6= l))．FLSA 推定の結果，同じ推定値をもつ区間によってデータは m = nF 個のクラスタに分割され，対応する説明変数のクラスタを C1 , · · · , Cm とする．ここで，FLSA 推定曲線に基づき，次の式で与えられるガウス型基底関数を構築することを考える．

}{(x − cj )2, (x < cj )exp −
2νh21jbj (x) =} (x = cj ) , j = 1, 2, · · · , m.
{ 1,2(x − cj ), (x > cj )
 exp − 2νh22j
指数関数 φ(x) = exp {−(x − µ)2 /(2h2 )} の 1 次導関数は(){}
x−µ(x − µ)20φ(x) = −exp −h22h22(4)(5)

であるため，x = µ − kh, µ + kh (k ≥ 1) における φ(x) の傾きが，それぞれ d1 ，d2 となる
ときの h の値 h1 , h2 は，( 2)( 2)
kkkk
h1 =
exp −, h2 =exp −|d1 |2|d2 |2(6)となる．従って，(6) 式より，x = µj ± kij hij における φ(x) の傾きが，FLSA 推定曲線の傾き dˆij と等しくなるような hij の値が次の式で与えられる．

kijĥij =exp|dˆij |()−kij2.
2(7)FF
ただし，xFj = {x1 j , · · · , xnjj } とし，dˆij は次のように定義する．

β̂F − β̂Fj−1 ˆβ̂F − β̂Fjdˆ1j = Fjj, d2j = Fj+1.
Fj−1F
xnj − xnj−1
x1 j+1 − x1 j(8)また，nj を j 番目のクラスタ Fj に含まれるデータ数とし，ガウス型基底関数の中心パラメータ cj を次の式で定義する．

ĉj =1 ∑xα .
nj x ∈Cα(9)j(4) 式の cj と hij をそれぞれ ĉj と ĥij で置き換えることによって, 以下の m 個の基底関数が得られる．

{}
2(x−ĉ)jexp −, (x < ĉj )2ν ĥ21jbj (x; ĉj , ĥ2j ) ={ 1,} (x = ĉj ) , j = 1, 2, · · · , m.
2(x − ĉj ), (x > ĉj ) exp −2ν ĥ22j(10)(10) 式の基底関数に基づく非線形回帰モデルは
yα = wT b(xα ) + α , α = 1, · · · , n,(11)で表され，誤差項が独立に平均 0，分散 σ 2 の正規分布に従うと仮定すると，密度関数は以下の式で与えられる．

]{yα − wT b(xα )}2f (yα |w, σ ) = √, α = 1, · · · , n.
exp −2σ 22πσ 221[3(12)

4

正則化最尤法
複雑な現象を近似するために非線形モデルを当てはめるとき，モデルのパラメータを最尤法によって推定すると，しばしばデータに過度に依存したモデルが推定される．そこ
で，u(x; w) の滑らかさが失われるにつれて増加する罰則項 (正則化項) を対数尤度関数に付け加えた
2lλ (w, σ ) =n∑log f (yα |w, σ 2 ) −

α=1nγ
H(w)
2(13)の最大化，すなわち正則化法によって推定する．ここで，γ (> 0) は正則化パラメータと呼ばれ，モデルの当てはまりのよさと曲線の滑らかさを制御する調整パラメータである．
また，H(w) は正則化項と呼ばれ，w の 2 次形式 w0 Kw などが用いられる．ただし K は既知の m + 1 次の非負定値行列である．このモデルのパラメータ w と σ 2 の正則化最尤推
定量は次で与えられる．

ŵ = (B T B + nγ σ̂ 2 K)−1 B T y,σ̂ 2 =1(y − B ŵ)T (y − B ŵ) .
n

ただし，B= (b(x1 )T , · · · , b(xn )T )T , y = (y1 , · · · , yn )T である．

5

モデル評価基準推定したモデルは，正則化パラメータ γ ，ガウス基底関数の幅調整パラメータ ν に依存しており，適切な値を決定することが重要となる．この問題に対し，モデル評価基準を用いて客観的に選択を行う．
統計的汎関数によって定義される一般の推定量の評価基準として，一般化情報量規準

GIC (Konishi and Kitagawa, 1996, 2008) がある．Konishi and Kitagawa (1996) の考えを，前述の正則化最尤法によって推定された基底関数展開に基づく非線形モデルに適用す
ると，回帰モデルにおける GIC はGIC = n(log 2π + 1) + n log σ̂ 2 + 2tr{R−1 Q}

4(14)で与えられる．ただし，

1 TT
21  B B + nγ σ̂ K σ̂ 2 B Λ1n R = ,n1 Tnσ̂ 21ΛBσ̂ 2 n2σ̂ 2
1 T 31 T
1 T 2T1  σ̂ 2 B Λ B − γK ŵ1n ΛB 2σ̂ 4 B Λ 1n − 2σ̂ 2 B Λ1n Q =
 1
1 T
1 T 4nT 3nσ̂ 21
Λ
B−1
ΛB1
Λ
1−n2σ̂ 4 n2σ̂ 2 n4σ̂ 6 n4σ̂ 2(15)(16)であり，1n は要素が全て 1 の n 次元ベクトル，Λ はΛ = diag(y1 − ŵT b(x1 ), · · · , yn − ŵT b(xn ))である．GIC の値が最小となる γ と ν の値を適切な値とする．

6

適用例
データは，yα = u(xα ) + v(xα ) + α から発生させ，真の曲線には次の関数を仮定した．
{u(x) = exp (−x)(x − 0.5)3 ,(a)v(x) = −1.5 · I(0.50 ≤ x ≤ 0.51),{u(x) = sin (2πx),(b)v(x) = −2 · I(0.25 ≤ x ≤ 0.26) + 2 · I(0.75 ≤ x ≤ 0.76),{u(x) = exp (−x) sin (5π exp (−x)),(c)v(x) = 0.

ただし， I(·) は定義関数とする．データ点 xα は区間 [0, 1] を等間隔に分割した点であり，誤差項 α に対しては平均 0，分散 τ 2 の正規分布を仮定し，τ = 0.2Ru と設定した．ただし，Ru は定義域 x ∈ [0, 1] に対する u(x) の値域の幅を表す．本節の数値実験では，λ1 = 0 と固定し，λ2 の適切な値を交差検証法によって選択した．また，正則化パラメータ γ と調整
パラメータ ν は，情報量規準 GIC で選択した．(13) 式の正則化項として H(w) = wT w を用いた．データ数 n = 100 に対して，FLSA に基づく提案手法 (FLSA-G) と，k-means 法によって構成された基底関数に基づく非線形回帰モデル (Km-G) と，FLSA 推定 (FLSA)の３つの手法を比較する実験を行った．
図 1 は真の曲線と推定曲線を比較したものである．表 1 は 50 回の数値実験に対する結果を示す．表中の mean m ，mean γ ，mean ν ，AMSE はそれぞれ，基底関数の個数 m の平均値，選択された正則化パラメータ γ の平均値，選択された調整パラメータ ν の平均，∑n
平均 2 乗誤差 MSE = α=1 {(u(xα ) + v(xα )) − ŷα }2 /n の平均値を表す．

5

1.5

0.10

11.0

0.05

0.50.00

y0.0

y0-0.20-1.0-0.15-1-0.5-0.10

y-0.0



0.0

0.2

0.4

0.6

0.8

1.0

0.00.2

x0.4

0.6

0.8

1.00.0

x0.2

0.4

0.6

0.8

1.0x

図 1: 真の曲線（点線）と，FLSA-G (実線), Km-G (鎖線)，FLSA (破線) の比較. 左上:(a). 右上: (b). 下: (c).

表 1: 数値実験結果．

function
(a)

(b)

(c)

FLSA-GKm-G
FLSA
FLSA-GKm-G
FLSA
FLSA-GKm-G
FLSA

mean m
8.20
6.72
−−
14.00
6.72
−−
13.20
7.54
−−

mean γ6.48 × 10−4
7.65 × 10−2−−
3.00 × 10−4
1.86 × 10−2−−
6.48 × 10−4
1.00 × 10−2−−

6mean ν
4.90
34.16
−−
402
2.22
−−
62.00
20.86
−−

AMSE
1.48 × 10−4
4.89 × 10−4
2.41 × 10−4
3.84 × 10−2
1.77 × 10−1
7.26 × 10−2
9.20 × 10−3
1.18 × 10−2
3.88 × 10−2

参考文献
[1] Ando, T., Konishi, S. and Imoto, S. (2008). Nonlinear regression modeling via regularized radial basis function networks. Journal of Statistical Planning and Inference
138, 3616–3633.
[2] Bishop, C. M. (1995). Neural Networks for Pattern Recognition. Oxford UniversityPress.
[3] de Boor, C. (2001). A Practical Guide to Splines. Springer.
[4] Eilers, P. and Marx, B. (1996). Flexible smoothing with B-splines and penalties (withdiscussion). Statist. Sci., 11, 89–121.
[5] Friedman, J., Hastie, T., Höfling, H. and Tibshirani, R. (2007). Pathwise coordinate
optimization. Ann. Appl. Statist., 1, 302–332.
[6] Hastie, T., Tibshirani, R. and Friedman, J. (2009). The Elements of Statistical Learning. 2nd ed. New York: Springer.
[7] Imoto, S. and Konishi, S. (2003). Selection of smoothing parameter in B-spline nonparametric regression models using information criteria. Annals of the Institute ofStatistical Mathematics 55, 671–687.
[8] Kawano, S. and Konishi, S. (2007). Nonlinear regression modeling via regularized
Gaussian basis functions. Bull. Inform. Cybern., 39, 83–96.
[9] Konishi, S., and Kitagawa, G. (1996). Generalised information criteria in modelselection. Biometrika 83, 875–890.
[10] Konishi, S. and Kitagawa, G. (2008). Information Criteria and Statistical Modeling.
Springer.
[11] Moody, J. and Darken, C. J. (1989). Fast learning in networks of locally-turnedprocessing units. Neural Computation. 1, 281–294.
[12] Tibshirani, R., Saunders, M., Rosset, S., Zhu, J. and Knight, K. (2005). Sparsity andsmoothness via the fused lasso. J. Roy. Statist. Soc. Ser. B, 67, 91–108.
7