これに関連して Lars の計算では以下を解く事が必要になります。
これらを XA に説明変数を加える、即ち列を加えながら解く事になります。
Eflon さんの論文では上を
と正規化しており、それを cholesky 分解を用いて解くものもあるようです(Zou and Hastie, 2005)。
効率を考えず、メモリも気にしなくていいなら正規化せず、例えばQR分解で解くことも可能ではないでしょうか？
ただ XA に一列加えて毎回一からQR分解を求めるのは効率が悪すぎでしょう・・・
そこで qrinsert、qrdelete を使う事が出来ます。
これらは Givens 回転を利用して前に求めた Q、R を更新して新たなQRを求めるものになります。
Givens回転は、以下の行列による座標回転を表します。
これにより行列の任意の要素を 0 にする事が出来ます。
上の G を左から掛けると行列の ( i, j ) 行がθだけ回転されます。
たとえば、行列
に Givens回転行列
を左から掛けると
となります。
ここで θを
となるように求めれば
と、３行１列を 0 に出来ます。
これを使い、XA に列を挿入する前に求めた Q、Rから挿入後の Q1、R1 を以下の様にして求める事が出来ます。
今、列の数が n である XA についてそのQR分解が
で、且つ XA の j 列目に x を挿入した行列のQR分解を
式（１）
式（２）
とします。
ここで R の j 列に Q' * x を挿入した
を考えると、これは XA1 の j 列の操作を飛ばして（j 列の対角成分以下を 0 にする Householder変換をしないで）求めた QR 分解と同じなので、この R0 でも
Q * R0 = XA1
が満たされます。
octave で試してみましょう。
XA と x を乱数で作ります。
octave:> XA = rand(4,3);
octave:> x = rand(4,1);
QR分解を求めます。
octave:> [Q,R]=qr(XA)
Q =
0.645464 -0.453811 -0.090253 -0.607689
0.069102 0.562767 0.690336 -0.449394
0.560067 0.655159 -0.475430 0.176232
0.514709 -0.219352 0.537827 0.630636
R =
0.97974 0.86334 0.50279
0.00000 0.77404 0.20891
0.00000 0.00000 0.21971
0.00000 0.00000 0.00000
XA1 として、 x を XA の 3列目に挿入したものを考えます。
octave:>XA1=[XA(:,1:2) x XA(:,3)]
XA1 =
0.632387 0.205988 0.457043 0.209893
0.067702 0.495263 0.117729 0.303990
0.548720 0.990648 0.248869 0.314006
0.504281 0.274582 0.569360 0.331131
R0 として R の3列目に Q'*x を挿入したものを計算します。
octave:> R0=[R(:,1:2) Q'*x R(:,3)]
R0 =
0.97974 0.86334 0.73558 0.50279
0.00000 0.77404 -0.10300 0.20891
0.00000 0.00000 0.22792 0.21971
0.00000 0.00000 0.07227 0.00000
octave:> Q*R0
ans =
0.632387 0.205988 0.457043 0.209893
0.067702 0.495263 0.117729 0.303990
0.548720 0.990648 0.248869 0.314006
0.504281 0.274582 0.569360 0.331131
XA1 と等しくなりました。
但し上の R0 では3列目以降上三角になっていません。
そこでこれを Givens回転で上三角にしてゆきます。
この為に R0 の ４行３列要素を 0 にするよう R0 の ３、４行を Givens回転します。
４行３列要素を 0 にする ３、４行の回転行列は
octave:> g=givens( R0(3,3), R0(4,3) );
３、４行のみ回転させる Givens回転行列は
octave:> G=[1 zeros(1,3); 0 1 zeros(1,2); zeros(2,2) g]
G =
1.00000 0.00000 0.00000 0.00000
0.00000 1.00000 0.00000 0.00000
0.00000 0.00000 0.95323 0.30226
0.00000 0.00000 -0.30226 0.95323
octave:> G*R0
ans =
0.97974 0.86334 0.73558 0.50279
0.00000 0.77404 -0.10300 0.20891
0.00000 0.00000 0.23910 0.20944
0.00000 0.00000 0.00000 -0.06641
これで R0 は対角化されました。
但し R に G を作用させたので、それを打ち消す操作を Q にも施さなければなりません。
G' * G = I
なので
XA1 = Q * R0 = ( Q * G' ) * ( G * R0 )
が成り立ちます。以上から
Q1 = Q * G'
R1 = G * R0
と、式（２）を満たす QR 分解が求められました。
アクティブセットから説明変数を外す、即ちXA から列を除く場合 (qrdelete) も R から対応する列を除いたものが上三角になるようGivens回転を施す事で同様に求められます。
このように一から QR 分解を求めずとも、前に求められた QR を元に、より少ない計算で新たな QR を求める事が出来ます。
この辺りの計算は octave の qrinsert、 qrdelete で行う事が出来ます。
octave:> [Q1, R1] = qrinsert(Q, R, 3, x)
Q1 =
0.645464 -0.453811 -0.269709 -0.551986
0.069102 0.562767 0.522214 -0.637033
0.560067 0.655159 -0.399926 0.311690
0.514709 -0.219352 0.703284 0.438578
R1 =
0.97974 0.86334 0.73558 0.50279
0.00000 0.77404 -0.10300 0.20891
0.00000 0.00000 0.23910 0.20944
0.00000 0.00000 0.00000 -0.06641
octave:> [Q2, R2] = qrdelete(Q, R, 3)
Q2 =
0.645464 -0.453811 -0.090253 -0.607689
0.069102 0.562767 0.690336 -0.449394
0.560067 0.655159 -0.475430 0.176232
0.514709 -0.219352 0.537827 0.630636
R2 =
0.97974 0.86334
0.00000 0.77404
0.00000 0.00000
0.00000 0.00000