#1. 要約
今回は重回帰分析のlassoペナルティを課した際の係数推定の記事です．特に係数ベクトルの要素ごとに１つずつ更新していく座標降下法(Coordinate Descent)について解説しています．実装に用いた言語はRです．本記事は，座標降下法の更新式の導出のみにフォーカスし，正則化法としてのlassoの数理などは扱っていません．
#2. lassoって？
lasso(Least Absolute Shrinkage and Selection Operator)はTibshirani(1996)で開発された線形回帰モデルの推定法で，以下の制約付き最小化問題を解くことを考えます．
\begin{eqnarray}
\underset{\beta}{min}\ &&\frac{1}{2n}||y-X\beta||^2_2=\frac{1}{2n}\sum_{i=1}^N(y_i-\sum_{j=1}^p\beta_jx_{ij})^2\\
&&s.t.\ ||\beta||_1\le s
\end{eqnarray}
ここで$||\cdot||^2_2$は$l_2$ノルムの二乗，$||\cdot||_1$は$l_1$ノルムを表しています．また，$X$は列中心化されており，正規化されているとします．適切なパラメータ$s$の下で以上の最小化問題を解くことによって，$\beta$のいくつかの要素をぴったり0として推定することができます．これはスパース推定と呼ばれる手法の最も基本的なもので，0となった部分に対応する独立変数は従属変数に寄与しないと解釈することができます．今回は座標降下法(Coordinate Descent)によって係数の推定を行います．座標降下法は，「スパース推定法による統計モデリング」では
$p$個あるパラメータのうち(p-1)個の要素の値を固定し，残る1個のパラメータのみを更新するという手順をすべての要素で繰り返すことにより更新値を得るアルゴリズムである．
(川野他「スパース推定法による統計モデリング」)
と説明されています．__ポイントは元の問題をパラメータ毎の最適化問題にを分割して，それぞれで最適化する__という点です．そのため，問題が分割できるかどうかをまず考えなければなりません．lassoの問題では以下のラグランジュ関数を最小化します．
\begin{eqnarray}
L_a(\beta) &=& \frac{1}{2n}||y-X\beta||^2_2 + a||\beta||_1\\
&=& \frac{1}{2n}\sum_{i=1}^N(y_i-\sum_{j=1}^p\beta_jx_{ij})^2+a\sum_{j=1}^p|\beta_j|
\end{eqnarray}
これをいつものように，$\beta$に関して偏微分して...といきたいところですが，$||\beta||_1$が微分できない点を持っています．そこで，その劣勾配(Subgradient)をまとめて
\begin{eqnarray}
d_j \in \left\{
\begin{array}{l}
{-1},\ \ (\beta_j<0)\\
[-1,1],\ \ (\beta_j=0)\\
{1},\ \ (\beta_j>0)
\end{array}\right.
\end{eqnarray}
と書くことにします．劣微分は微分の概念を拡張したものだと思っておいてください．ポイントは$\beta_j=0$の時に，$[-1,1]$のように範囲として値をとることです．この$d_j$の記法を用いて，ラグランジュ関数を微分して0とおくと，
\begin{eqnarray}
\frac{\partial}{\partial \beta_j}L_a(\beta)&=&-\frac{1}{n}\sum_{i=1}^Nx_{ij}(y_i-\sum_{k=1}^p\beta_kx_{ik})+ad_j=0\\
&\Leftrightarrow& \frac{1}{n}\sum_{i=1}^Nx_{ij}(y_i-\sum_{k\ne j}^p\beta_kx_{ik}-\beta_jx_{ij})=ad_j\\
&\Leftrightarrow& \beta_j\frac{1}{n}\sum_{i=1}^Nx_{ij}^2=\frac{1}{n}\sum_{i=1}^Nx_{ij}(y_i-\sum_{k\ne j}^p\beta_kx_{ik})-ad_j\\
&\Leftrightarrow& \beta_j=\frac{1}{n}\sum_{i=1}^Nx_{ij}(y_i-\sum_{k\ne j}^p\beta_kx_{ik})-ad_j\tag{1}\\
&&(\because X\mbox{は正規化されており，}n^{-1}\sum_{i=1}^Nx_{ij}^2=1)
\end{eqnarray}
という$\beta_j$に関する推定方程式が得られます．次に$\beta_j$以外のパラメータを固定したものを，$\tilde{\beta_k}(k\neq j)$とおいて，第j変数に対応する項を除いた残差ベクトルを
\begin{eqnarray}
r^{(j)} = y - \sum_{k\ne j}x_{(k)}\tilde{\beta_k}
\end{eqnarray}
とおくと，(1)式は
\begin{eqnarray}
\beta_j=\frac{1}{n}\sum_{i=1}^Nx_{ij}r^{(j)}_i-ad_j
\end{eqnarray}
とスッキリした形で書き表せます．さらに，以下で定義される
\begin{eqnarray}
S(x,\lambda) &=& sign(x)(|x|-\lambda)_{+}
\end{eqnarray}
軟閾値作用素(Soft thresholding operator)を用いると，
\begin{eqnarray}
\beta_j = S(\frac{1}{n}r^{(j)T}x_{(j)},a)
\end{eqnarray}
と表されます．座標降下法では，$\beta$を要素が収束するまで更新することを繰り返します．今回のプログラムでは，更新された$\beta$とその前のステップの$\beta$の要素の差の中で最も大きいものが$10^{-3}$より小さくなった時に収束したと判断しています．
#3. プログラム
今回は，Rで実装しています(プロフィールにRと書いているのに，Rのスクリプト全然あげてないから&授業で書いたのを使いまわせるから)．探せばlassoのパッケージあるので，実データを解析する際は，そちらを使いましょう．解析に用いたデータは，https://web.stanford.edu/~hastie/StatLearnSparsity/data.htmlの米国犯罪データです．
soft.th = function(lambda,x){
return(sign(x)*pmax(abs(x)-lambda,0))
}
lasso = function(X, y, lambda=0){
X=as.matrix(X); X=scale(X); p=ncol(X); n=length(y); X.bar=array(dim=p);
for(j in 1:p){X.bar[j]=mean(X[,j]);X[,j]=X[,j]-X.bar[j];};
y.bar=mean(y); y=y-y.bar
eps=1; beta=array(0, dim=p); beta.old=array(0, dim=p)
while(eps>0.001){
for(j in 1:p){
r= y - X[,-j]%*%beta[-j]
beta[j]= soft.th(lambda,cov(r,X[,j]))
}
eps=max(abs(beta-beta.old)); beta.old=beta
}
beta.0= y.bar - X.bar%*%beta
return(list(beta=beta, beta.0=beta.0))
}
df=read.table("crime.txtへの適当なパス"); x=df[,3:7]; y=df[,1];
lasso(x,y,40)
Rは統計処理に特化した言語で，今回のようにスクラッチで実装するよりも，便利な関数を使って実データ分析に用いるのが良いと思います．
#参考文献
川野 秀一・松井 秀俊・廣瀬 慧(2018)「スパース推定法による統計モデリング」
Friedman, Hastie, Tibshirani(2010). Regularization paths for generalized linear models via coordinate descent