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 ! ループの１回目と２回目はここをスキップ alm=alm*alf ! ループの３回目からは alm を alf を乗じる goto 10291 10311 if(m .ne. 1)goto 10321 ! ループの２回目はここをスキップ alm=big ! ループの１回目は alm = big(9.9d35) にする goto 10331 10321 continue alm=0.0 ! ループの２回目は 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になると思われます。
! ２番目のループ ! 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 ! ２番目のループここまで
上記の処理で変数を横断して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
ちょっと長くなってしまったので一度きります。 次回はループ③から始めます。
