glmnetをもう少し理解したい②

前回の記事では glmnet の中身を確認し、引数の family によって呼び出す関数を変えていることがわかりました。 今回はそのなかでも gaussian が指定された場合の関数である elnet を見ていきましょう。 なお前回の記事はこちらです。

ushi-goroshi.hatenablog.com

elnet の実装

それでは早速 elnet という関数を見ていきましょう。 ちなみにここでの elnet はコンソールで elnet と打っても表示されませんが、C ã‚„ Fortran で書かれたものではなくて単に glmnet からエクスポートされていない関数なので glmnet:::elnet で中身を見ることができます。

この関数はそれほど長くないのでいきなり内容の確認に入りますが、他の多くの関数同様に elnet でも最初はパラメータの受け取り・確認を行います。 下のブロックでは反復回数（ maxit ）、観測値の重み（ weights ）を受け取った後、 type.gaussian の指定内容によって ka というパラメータに格納する値を変えています。

function (x, is.sparse, ix, jx, y, weights, offset, type.gaussian = c("covariance", 
    "naive"), alpha, nobs, nvars, jd, vp, cl, ne, nx, nlam, flmin, 
    ulam, thresh, isd, intr, vnames, maxit) 
{
    # 1. パラメータの受け取り
    ### maxit
    maxit = as.integer(maxit)
    ### weights
    weights = as.double(weights)
    ### type.gaussian
    type.gaussian = match.arg(type.gaussian)
    ka = as.integer(switch(type.gaussian, covariance = 1, naive = 2, 
        ))

ka はさらに先の処理で elnetu と elnetn という２つのサブルーチンのどちらを呼ぶかを決めていますので、 type.gaussian の指定に合わせてサブルーチンを変更しているということですね。

以下では y および offset （存在する場合）を double に変換しています。 また y の重み付き平均を使って Null Deviance （残差逸脱度）を計算しています。

### y の storage.mode
storage.mode(y) = "double"
### offset
if (is.null(offset)) {
  is.offset = FALSE
}
else {
  storage.mode(offset) = "double"
  is.offset = TRUE
  y = y - offset
}
### 重み付き平均
ybar = weighted.mean(y, weights)
### Null Deviance（帰無モデルの残差逸脱度）
nulldev = sum(weights * (y - ybar)^2)
if (nulldev == 0) 
  stop("y is constant; gaussian glmnet fails at standardization step")

次のブロックで早速フィッティングに入ります。 is.sparse が指定されているか否かで spelnet と elnet のどちらが呼ばれるかが決まりますが、引数の違いとしては spelnet において x が as.double とされており、 ix と jx （いずれも疎行列において非ゼロの要素の座標を特定するための数値）が追加されています。

# 2. フィッティング
## 疎行列であるかで関数を変える
fit = if (is.sparse) 
  .Fortran("spelnet", ka, parm = alpha, nobs, nvars, x, 
           # 疎行列である場合、以下の ix, jx が引数として追加される
           # ix, jx は疎行列における非ゼロの要素の累積個数と行番号
           ix, jx, 
           y, weights, jd, vp, cl, ne, nx, nlam, flmin, 
           ulam, thresh, isd, intr, maxit, lmu = integer(1), 
           a0 = double(nlam), ca = double(nx * nlam), ia = integer(nx), 
           nin = integer(nlam), rsq = double(nlam), alm = double(nlam), 
           nlp = integer(1), jerr = integer(1), PACKAGE = "glmnet")
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 = 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 の個数
# なので ca は変数の数 * lambda の数

処理を抜けたあとは、エラーをチェックした上で必要なパラメータを取得します。

# 3. 後処理
## エラーチェック
if (fit$jerr != 0) {
  errmsg = jerr(fit$jerr, maxit, pmax = nx, family = "gaussian")
  if (errmsg$fatal) 
    stop(errmsg$msg, call. = FALSE)
  else warning(errmsg$msg, call. = FALSE)
}
## パラメータ（切片、回帰係数、自由度、次元、lambda）を取ってくる
outlist = getcoef(fit, nvars, nx, vnames)

## パラメータ（xxxxxxxxxx）を取ってきて outlist に結合する
dev = fit$rsq[seq(fit$lmu)]

outlist = c(outlist, list(dev.ratio = dev, nulldev = nulldev, 
                          npasses = fit$nlp, jerr = fit$jerr, offset = is.offset))
## elnet クラスを付与する
class(outlist) = "elnet"
outlist
}

それでは次に elnet の本体である elnet（ややこしいですね） の中身を見ていきましょう。

elnet（二度目）の実装

上記のフィッティングのセクションで elnet は .Fortran("elnet") として呼ばれていました。これまで glm ã‚„ GAM で見てきたときと同じように、 glmnet でもやはり fortran に行き着くようですね。

と言ってもここではまだ関数自体は大きくなく、下のように（コメント抜きで）30行程度で書かれています。

subroutine elnet(ka,parm,no,ni,x,y,w,jd,vp,cl,ne,nx,nlam,  flmin,u
                 *lam,thr,isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
implicit double precision(a-h,o-z)                                
double precision x(no,ni),y(no),w(no),vp(ni),ca(nx,nlam),cl(2,ni) 
double precision ulam(nlam),a0(nlam),rsq(nlam),alm(nlam)          
integer jd(*),ia(nx),nin(nlam)                                    
double precision, dimension (:), allocatable :: vq;               

! vp が 0.0 だった場合には jerr = 100000 として return してしまう
if(maxval(vp) .gt. 0.0)goto 10021                                 
jerr=10000                                                        
return                                                            
10021 continue                                                          
allocate(vq(1:ni),stat=jerr)

! ここでも jerr に 0 以外の数値が入っていたら return してしまう
if(jerr.ne.0) return                                              

! vp の値によって vq を生成
! デフォルトは 1
! ni は nvars で変数の数なので、 vq にはデフォルトでは変数の数が入る
! でもなんで sum(vq) なんだろ
vq=max(0d0,vp)                                                    
vq=vq*ni/sum(vq)

! elnetu か elnetn のどちらを呼ぶかは ka .ne. 1 であるかで判断している
! 1 でなければ elnetn 、 1 なら elnetu
if(ka .ne. 1)goto 10041                                           
call elnetu  (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,
              *isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
goto 10051                                                        
10041 continue                                                          
call elnetn (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,i
             *sd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
10051 continue                                                          
continue                                                          
deallocate(vq)                                                    
return                                                            
end                                                               

goto を多用していますね。。。 変数宣言以下で気になるところとしては、 vp が 0 だったときの挙動と、 elnetu を呼ぶところでしょうか。

vp は前回の記事で確認した通り、 glmnet のなかで vp = as.double(penalty.factor) として定義されています。 この penalty.factor はデフォルトでは 1 が入りますので基本的には goto 10021 で飛ばされてしまいます。 このセクションで引っかかるのは明示的に penalty.factor に 0 を指定した場合ですね。

! vp は各変数に対する罰則の重み（デフォルトは 1） が入ったベクトル
! vp = as.double(penalty.factor)
! jerr の数値で後続の処理で出力するエラーメッセージが決まる
if(maxval(vp) .gt. 0.0)goto 10021                                 
jerr=10000                                                        
return                                                            
10021 continue                                                          
allocate(vq(1:ni),stat=jerr)

では penalty.factor に 0 を指定した場合はどうなるかと言うと、 jerr に 10000 が入力されて return されます。 この jerr は先ほど確認した後処理において errmsg = jerr(fit$jerr, maxit, pmax = nx, family = "gaussian") としてエラーメッセージに変換されるのでした。 またこの jerr という関数は glmnet で定義されていますので、

> glmnet:::jerr
function (n, maxit, pmax, family) 
{
  if (n == 0) 
    list(n = 0, fatal = FALSE, msg = "")
  else {
    errlist = switch(family, gaussian = jerr.elnet(n, maxit, 
                                                   pmax), binomial = jerr.lognet(n, maxit, pmax), multinomial = jerr.lognet(n, 
                                                                                                                            maxit, pmax), poisson = jerr.fishnet(n, maxit, pmax), 
                     cox = jerr.coxnet(n, maxit, pmax), mrelnet = jerr.mrelnet(n, 
                                                                               maxit, pmax))
    names(errlist) = c("n", "fatal", "msg")
    errlist$msg = paste("from glmnet Fortran code (error code ", 
                        n, "); ", errlist$msg, sep = "")
    errlist
  }
}

として取り出せます。 関数をみてみると、 errlist は switch(family, ~) で更に異なる関数を呼び出し、その結果を格納しているようです。 そのため更に jerr.elnet を確認すると

> glmnet:::jerr.elnet
function (n, maxit, pmax) 
{
  if (n > 0) {
    if (n < 7777) 
      msg = "Memory allocation error; contact package maintainer"
    else if (n == 7777) 
      msg = "All used predictors have zero variance"
    else if (n == 10000) 
      msg = "All penalty factors are <= 0"
    else msg = "Unknown error"
    list(n = n, fatal = TRUE, msg = msg)
  }
  else if (n < 0) {
    if (n > -10000) 
      msg = paste("Convergence for ", -n, "th lambda value not reached after maxit=", 
                  maxit, " iterations; solutions for larger lambdas returned", 
                  sep = "")
    if (n < -10000) 
      msg = paste("Number of nonzero coefficients along the path exceeds pmax=", 
                  pmax, " at ", -n - 10000, "th lambda value; solutions for larger lambdas returned", 
                  sep = "")
    list(n = n, fatal = FALSE, msg = msg)
  }
}

else if (n == 10000) msg = "All penalty factors are <= 0" と、罰則項が 0 であることを教えてくれていますね。

さて続いて elnetu の呼びだしを確認すると、elnetu と elnetn のいずれを呼ぶかは ka で決まっています。 先ほど少し触れた通り、 ka は ka = as.integer(switch(type.gaussian, covariance = 1, naive = 2, )) で定義されています。 また type.gaussian は glmnet の引数であり、type.gaussian = ifelse(nvars < 500, "covariance", "naive") と定義されています。 変数の数が 500 未満であれば covarinace となり、 ka には 1 が引き渡されるので if(ka .ne. 1) には該当せず、したがって elnetu が呼ばれることになるようですね。

! elnetu か elnetn のどちらを呼ぶかは ka .ne. 1 であるかで判断している
! 1 でなければ elnetn 、 1 なら elnetu
! ka は elnet の第一引数
! ka = as.integer(switch(type.gaussian, covariance = 1, naive = 2, ))
! この covariance / naive は変数の数で決まる
! type.gaussian = ifelse(nvars < 500, "covariance", "naive")
if(ka .ne. 1)goto 10041                                           
call elnetu  (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,
              *isd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)
goto 10051                                                        
10041 continue                                                          
call elnetn (parm,no,ni,x,y,w,jd,vq,cl,ne,nx,nlam,flmin,ulam,thr,i
             *sd,intr,maxit,  lmu,a0,ca,ia,nin,rsq,alm,nlp,jerr)

次回はこの elnetu を見てみましょう。