glmnetをもう少し理解したい④
それでは前回の記事に続いてelnet1
の紹介です。過去の記事はこちらです。
ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com ushi-goroshi.hatenablog.com
elnet1
の実装
前回の記事で最後に触れた通り、elnet1
自体は 180 行程度とそれほど大きくはないサブルーチンなのですが、多数のループが込み入っています。
具体的には以下の通り 9 つのループ処理(fortran なので do 文)がネストした構造となっており、 しかもgoto
によって行き来しています(わかりやすいように R で書いてありますが、添字は統一してあります)。
# 1番目 for (m in 1:nlam) { # 2番目のループ for (j in 1:ni) { } # 3番目のループ for (k in 1:ni) { # 4番目のループ for (j in 1:ni) { } # 5番目のループ for (j in 1:ni) { } } # 6番目のループ for (l in 1:nin) { # 7番目のループ for (j in 1:nin) { } } # 8番目のループ for (j in 1:ni) { } # 9番目のループ for (j in 1:nin) { } }
前処理
まずはいつもの通り変数の定義ですが、それに加えて初期パラメータを取得するという処理が入ります。
subroutine elnet1(beta,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,ulam,th *r,maxit,xv, lmu,ao,ia,kin,rsqo,almo,nlp,jerr) implicit double precision(a-h,o-z) double precision vp(ni),g(ni),x(no,ni),ulam(nlam),ao(nx,nlam) double precision rsqo(nlam),almo(nlam),xv(ni) double precision cl(2,ni) integer ju(ni),ia(nx),kin(nlam) double precision, dimension (:), allocatable :: a,da integer, dimension (:), allocatable :: mm double precision, dimension (:,:), allocatable :: c allocate(c(1:ni,1:nx),stat=jerr) if(jerr.ne.0) return; ! 初期パラメータを取得 call get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx,itrace) ! a, mm, da を allocate allocate(a(1:ni),stat=jerr) ! a は説明変数の数の次元をもつベクトル if(jerr.ne.0) return allocate(mm(1:ni),stat=jerr) ! mm は説明変数の数の次元をもつベクトル if(jerr.ne.0) return allocate(da(1:ni),stat=jerr) if(jerr.ne.0) return
ここでget_int_parms
はそれほど大きくないので全体を見てみましょう。
以下のようなサブルーチンです:
subroutine get_int_parms(sml,eps,big,mnlam,rsqmax,pmin,exmx,itrace) implicit double precision(a-h,o-z) data sml0,eps0,big0,mnlam0,rsqmax0,pmin0,exmx0,itrace0 /1.0d-5,1.0d-6,9.9d35,5,0.999,1.0d-9,250.0,0/ sml=sml0 eps=eps0 big=big0 mnlam=mnlam0 rsqmax=rsqmax0 pmin=pmin0 exmx=exmx0 itrace=itrace0 return entry chg_fract_dev(arg) sml0=arg return entry chg_dev_max(arg) rsqmax0=arg return entry chg_min_flmin(arg) eps0=arg return entry chg_big(arg) big0=arg return entry chg_min_lambdas(irg) mnlam0=irg return entry chg_min_null_prob(arg) pmin0=arg return entry chg_max_exp(arg) exmx0=arg return entry chg_itrace(irg) itrace0=irg return end
上から3行目のdata
文は変数に初期値を与える fortran の記法のようで、data
に続いて宣言した変数に対して/
で挟んだ値を初期値として与えるようです。
そのためsml0
には 1.0d-5 が、eps0
には 1.0d-6 が入力されます。
ここで d は倍精度の指数表記を表します。13行目のentry
以降は各変数について特定の値を指定するためのもののようです(entry
の使い方がよくわからない…)。
続けていくつか変数に値を代入します。
まずはbta
ですが、代入しているbeta
は元々parm
として渡されたもので、これはelnet.r
でparm = alpha
として渡していたものでした。さらにこのalpha
はglmnet.r
で定義されたもので、L1 と L2 それぞれに対する罰則の配分を決めるパラメータです:
(なぜか Tex が表示されないのでひとまず)
(1 − \alpha)/2 ||\beta||^2_2 + \alpha||\beta||_1
bta=beta
このbta
を 1 から減じたものをomb
としますが、このomb
はすぐ下で定義されるalm
との乗算でdem
を定義する(つまりdem = alm * obm
)ためだけに使われています。
さらにalm
はループの中で更新されながら最終的にはbta
との乗算によってab
となり、回帰係数の縮小に使われることになります。
またその次のalf
はalm
の更新に使われますので、これらの変数がループの中で更新されつつ回帰係数の縮小に利用されるということになります(他にもあります)。
omb=1.0-bta alm=0.0 alf=1.0
以下のブロックではeqs
とalf
を定義しますが、flmin
が 1.0以上であればスキップされるようです。
このflmin
というのはglmnet.r
において罰則lambda
が指定されていれば 1
が、されていない時にはlambda.min.ratio
が入力される変数でした。
lambda.min.ratio
はデフォルトではlambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04)
となっていますので 1 よりは小さい値が入りそうです。
したがって以下のブロックは「lambda
が指定されていないときはalf
を定義しよう」という処理になっています(eqs
はここしか出てきません)。
その場合、eps
とflmin
(=1)の大きい方を新たにeqs
と定義しますが、このeps
はget_int_parms
でeps0
(1.0d-6 という小さい数)を受け取っていました。
一方lambda.min.ratio
は先ほど述べたようにデフォルトではlambda.min.ratio = ifelse(nobs < nvars, 0.01, 1e-04)
となっていますので、もう少し大きい値となりそうです。
したがってeqs
は 0.01 or 1e-04 、alf
はその1/(nlam-1)
乗となるようです。
if(flmin .ge. 1.0)goto 10271 eqs=max(eps,flmin) alf=eqs**(1.0/(nlam-1)) ! alf を eqs の (1/(nlam-1)) で定義する
flmin
が 1 以上である(lambda
が指定されている)場合は上記をスキップしてこちらにきます。rsq
はそのまま残差平方和ですね。
続くa
はelnet1
の中で重要な役割を担っているのでじっくりと見ていきましょう。
実はこのa
は(縮小された)回帰係数を格納する変数です。
10271 continue ! パラメータの初期化 rsq=0.0 ! 残差平方和 a=0.0
このa
がどうなるのか、フライングして先の処理を見てみましょう。
elnet1
の 70 行目前後に以下の処理があります:
ak=a(k) u=g(k)+ak*xv(k) v=abs(u)-vp(k)*ab a(k)=0.0 if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*dem))) if(a(k).eq.ak)goto 10371
ak
という変数にa
の k 番目の値を渡しておき、u
とv
を定義し、a
の k 番目の値を 0 に更新した上で色んな値を参照しながら再度更新しています(このu
やv
は後で確認します)。
最終的にa
は以下のようにao
という変数に代入されます(154 行目):
if(nin.gt.0) ao(1:nin,m)=a(ia(1:nin))
このao
ですが、elnetu
の中でelnet1
を呼び出すときにはca
という引数として渡されています。
! elnet1 で受け取る変数 ! lmu の次に ao がある subroutine elnet1(beta,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,ulam,th *r,maxit,xv, lmu,ao,ia,kin,rsqo,almo,nlp,jerr) ! elnetu で elnet1 を call するときの引数 ! こちらは lmu の次に ca がある call elnet1(parm,ni,ju,vp,cl,g,no,ne,nx,x,nlam,flmin,vlam,thr,maxi,xv, lmu,ca,ia,nin,rsq,alm,nlp,jerr)
このca
はelnet.r
の中で.Fortran("elnet", ...)
と call される際に定義される変数でした:
else .Fortran("elnet", ka, parm = alpha, nobs, nvars, as.double(x), y, weights, jd, vp, cl, ne, nx, nlam, flmin, ulam, thresh, isd, intr, maxit, lmu = integer(1), a0 = double(nlam), # ここで ca が定義されている ca = double(nx * nlam), ia = integer(nx), nin = integer(nlam), rsq = double(nlam), alm = double(nlam), nlp = integer(1), jerr = integer(1), PACKAGE = "glmnet")
ここでnx
は説明変数の数、nlam
は罰則lambda
の数なので、説明変数の数 × lambda
の数のベクトルを定義しています(そしてそれがelnet1
の中でao
として評価・格納される)。
このca
はelnet.r
の後続の処理において以下の箇所で抽出されます:
outlist = getcoef(fit, nvars, nx, vnames)
ここでglmnet:::getcoef
は以下の通りで、fit
として返ってきたオブジェクトのca
そのものをbeta
に格納しています(ninmax
が 0 の場合は 0 のベクトルが返る)。
# glmnet:::getcoef function (fit, nvars, nx, vnames) { # ここまで省略 nin = fit$nin[seq(lmu)] ninmax = max(nin) # ここまで省略 if (ninmax > 0) { # ここで ca を抽出している ca = matrix(fit$ca[seq(nx * lmu)], nx, lmu)[seq(ninmax), , drop = FALSE] df = apply(abs(ca) > 0, 2, sum) ja = fit$ia[seq(ninmax)] oja = order(ja) ja = rep(ja[oja], lmu) ia = cumsum(c(1, rep(ninmax, lmu))) # beta に格納する beta = drop0(new("dgCMatrix", Dim = dd, Dimnames = list(vnames, stepnames), x = as.vector(ca[oja, ]), p = as.integer(ia - 1), i = as.integer(ja - 1))) } else { beta = zeromat(nvars, lmu, vnames, stepnames) df = rep(0, lmu) } # ここも省略 list(a0 = a0, beta = beta, df = df, dim = dd, lambda = lam) }
これにいくつかの情報を追加したものがglmnet
の返り値です。elnet1
において評価されたa
がao
に格納され、elnet
にca
として渡され、elnet.r
でbeta
に抽出・格納される流れが伝わりましたでしょうか。
重要な変数を説明したところなので、以下ブロックで初期化している変数の詳細は出てきたときに説明するとして、さっさと次に進んでしまいましょう。
mm=0 nlp=0 nin=nlp iz=0 mnl=min(mnlam,nlam)
ループ①(alm
の更新)
上記までで必要な変数の初期化が完了したので、以下よりループに入ります。
一番外側のループはlambda
の個数(nlam
)に対して実行されますが、nlam
のデフォルトは 100 となっています(glmnet.r
)。
以下ではおおよそalm
を更新する処理を行うのですが、lambda
の指定の有無や、ループの回数によってalm
に入力する値を変えています。
まずはlambda
の指定の有無で処理を分けます。以下のまとまりはflmin
が 1.0 より小さい場合にスキップされますが、先ほど述べたように、flmin
はglmnet.r
においてlambda
の指定がない場合に相当します。
lambda
の指定がある場合にはalm = ulam(m)
としてalm
を更新した上で、10291 までスキップするのですが、この 10291 は 2 番目のループの中にありますので、少し大きめのスキップとなるようです。
なおulam
はlambda
が指定されている場合、lambda
の降順になっているため、ループの 1 回目であればlambda
の最大値が入ります。
do 10281 m=1,nlam ! nlambda なので lambda の個数だけループ if(itrace.ne.0) call setpb(m-1) ! プログレスバー if(flmin .lt. 1.0)goto 10301 alm=ulam(m) ! flmin が 1.0 以上の場合は alm = ulam(m) とする goto 10291
lambda
の指定がなければ以下の処理に入るのですが、ここではループの回数によってalm
に入力する値を変えています。
具体的には、ループの 1 回目にはbig
(9.9d35)という極端に大きな値を入力し、 2 回目には 0.0 を、3 回目以降は 元の値にalf
を乗じたものを入力します。
10301 if(m .le. 2)goto 10311 ! ループの1回目と2回目はここをスキップ alm=alm*alf ! ループの3回目からは alm を alf を乗じる goto 10291 10311 if(m .ne. 1)goto 10321 ! ループの2回目はここをスキップ alm=big ! ループの1回目は alm = big(9.9d35) にする goto 10331 10321 continue alm=0.0 ! ループの2回目は alm を いったん 0 にする
このalf
は先ほど説明した通りeqs^(1.0/(nlam-1))
として定義されますが、eqs
が 0.01 or 1e-4 とすると、nlambda
を 10 とした場合には以下のような数値になります:
0.01^(1/(10-1)) # [1] 0.5994843 1e-4^(1.0/(10-1)) # [1] 0.3593814
つまりalm
はだんだん絶対値が小さくなるわけですね。
ループ②(罰則の定義)
続いて 2 番目のループに入ります…と言いつつ 2 番目のループは一瞬で終わります。
先ほど更新したalm
について変数ごとの内積と比較し、大きい方を採用します。
したがってここでは各変数に対するループとなります。
まずju
とvp
ですが、ju
は前回記事で確認した通り、chkvars
によって各変数列の内容が全く同じでないかを確認したものでした。
ある変数列の中身が全く同じであれば 0 であったため、ここで次の変数にスキップされます。
次にvp
ですが、これは 1 回目の記事で確認した通りglmnet.r
において各変数に対する罰則の重み(デフォルトは 1) が入ったベクトルとして定義されたものです(vp = as.double(penalty.factor)
)。
罰則をかけない場合は 0 となり、スキップされるようです。
変数にバラつきがあり、罰則を検討する場合にはここで再度alm
を更新します。
ここで出てくるg
はstandard
の中でy
とx
の内積(共分散)を格納したものとして定義されたものでした。
それを罰則の大きさで除しているため、penalty.factor
を小さくする(分母が小さくなる)と共分散が大きくなり変数として残りやすい、というロジックになっているようですね。
ちなみにループ①の 1 回目のループはalm
に 9.9d35 という数値が入るので必ずこの数値が採用されると思います。またループ 2 回目は今度はalm
が 0.0 になるため、今度は必ず変数の共分散側の数値がalm
になると思われます。
! 2番目のループ ! alm の更新 do 10341 j=1,ni ! ni は変数の数 if(ju(j).eq.0) goto 10341 if(vp(j).le.0.0) goto 10341 alm=max(alm,abs(g(j))/vp(j)) 10341 continue ! 2番目のループここまで
上記の処理で変数を横断してalm
を更新したのち、以下でさらにalm
を更新します。
ここではbta
(alpha
; L1 と L2 への重みの配分パラメータ)と 0.001 の max で alm
を除し、alf
を乗じています。
一応ここで式を確認しておくと以下のようになります:
一体これは何をやっているんでしょうか…。
continue alm=alf*alm/max(bta,1.0d-3)
続いていくつかの変数を更新します。
dem
はalm * omb
として定義されますが、ここでomb
は (1-bta
)でした。 またab
はalm
にbta
を乗じたものですので、これらはそれぞれ「lambda
×(1-alpha)
」および「lambda
×alpha
」ということになり、dem
とab
が実質的な罰則の大きさを表すことになりそうですね。
10331 continue 10291 continue dem=alm*omb ! dem = alm * (1-bta) ab=alm*bta ! ab = alm * bta
これらがどのように使われているか少し先を見てみましょう。
! ab u=g(k)+ak*xv(k) ! L69(ループ③の中)、L119(ループ⑥の中) v=abs(u)-vp(k)*ab ! L70、L120(ともに上に同じ) ! dem a(k)=0.0 ! L71、L121 if(v.gt.0.0) a(k)=max(cl(1,k),min(cl(2,k),sign(v,u)/(xv(k)+vp(k)*dem))) ! L72、L122
両方ともvp
に乗じており、ab
はabs(u)
からの減算、dem
はxv(k)
との加算の後にsign(v,u)
と除算し、cl
との max/min を取っています。
vp
は罰則の重みを定義したものでしたので、alpha
とlambda
で決まる罰則の大きさをそのまま使うか弱くするかを決めています。
dem
の方は演算の結果をa
に格納していますが、前述の通りa
は回帰係数を保存する変数でしたので、sign(v,u)/(xv(k)+vp(k)*dem)
がcl(1,k)
よりも大きければa
、すなわち回帰係数が更新されるということになりますね。
またこの演算が実行されるかの基準としてv
が使われており、このv
を計算するためにab
が使われている、ということのようです。
じゃあこのu
とかv
って何なの?という話なのですが、これは次のループの話なので少しお待ちください。
残る変数のうちrsq0
は残差平方和ですね。またjz
はiz
と組み合わせて使われていますが、この条件分岐がちょっと理解出来なかったのでスキップします。
一応、iz
はループ①の途中(ループ③が終了した時点)で 1 になるため、iz * jz
が 0 になるのはほぼjz
が 0 の時に限ると言えそうです。
nlp
は iteration のカウンターとして使われており、dlx
は回帰係数の更新前後の差分を見ています。
どちらもループを抜けるための基準として使われています。
rsq0=rsq jz=1 continue 10351 continue if(iz*jz.ne.0) goto 10360 ! iz = 0, jz = 1 nlp=nlp+1 dlx=0.0
ちょっと長くなってしまったので一度きります。 次回はループ③から始めます。