R の線形回帰 lm 関数が完全な多重共線性のある入力データをどう処理しているかを辿る
畳屋民也やること
R で default で利用可能な線形回帰の関数 lm は、完全な多重共線性のあるデータ入力に対して、特定の説明変数を除去して実行結果を返す（singular.ok=FALSE を指定すると error を吐いて止まる）。
このとき、どのような基準で残す変数を決定しているかをコードをたどりながら調べる。
畳屋民也lm のコードはどこで読めるか？
そもそも lm は stats というパッケージの中に含まれている。
このうち、lm の実装は以下：
回帰係数を最小二乗法により推定する際に、QR 分解を用いているらしい。
そして、QR 分解については C 言語で実装されていそうである：
→ ここからさらに FORTRAN で実装された関数を呼び出していた。
畳屋民也QR 分解と最小二乗推定 (OLS)
とりあえず見つけた資料を貼る。
QR 分解については、
畳屋民也愚直に lm のコードを追ってみる
主に回帰係数の推定部分の処理がどこで行われているかを見て見た。
[https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/lm.R]
lm
変数 z に fit の結果が入ってくる。
推定部分は lm.fit() or lm.wfit() を呼び出している（重み weights を指定しているかで分岐）。
lm.fit
QR 分解アルゴリズムは C_Cdqrls を呼び出している。
Cdqrls
呼び出している関数 C_Cdqrls は、 C 言語で書かれた こちら で定義されている関数 Cdqrlsと考えられる。
なお、コメントに書かれているように、FORTRAN で書かれた関数の wrapper らしい。
最終的には下記で dqrls という FORTRAN の関数を呼び出している。
dqrls
[https://github.com/SurajGupta/r-source/blob/master/src/appl/dqrls.f]
この中から dqrdc2 と dqrsl という別の関数を呼び出している：
dqrdc2
[https://github.com/SurajGupta/r-source/blob/master/src/appl/dqrdc2.f]
Householder Transformation という方法で QR 分解をおこなっているらしい。
dqrsl
[https://github.com/SurajGupta/r-source/blob/master/src/appl/dqrsl.f]
dqrdc2 の出力結果を用いて最小二乗法による回帰係数を算出しているものと思われる。
その他
FORTRAN を真面目に読むのであれば、下記が参考になりそう。
畳屋民也入力した説明変数に完全な多重共線性がある場合、どこで変数選択をしているか？
計画行列を Householder Transformation により QR分解する際に、説明変数間で完全な多重共線性が見つかれば除外している。
コードでいうと、以下の中で処理が行われている。
[https://github.com/SurajGupta/r-source/blob/master/src/appl/dqrdc2.f]
より詳細な話
QR分解を行う際、各ステップごとに Pivoting といって列の入れ替えを行なっている。
これにより、l2ノルムが大きい列（≒ それまでに求めた基底ベクトルに射影した後の残成分が大きい説明変数）から順に処理が行われる。
完全な多重共線性がある（計画行列Xが full-rank でない）場合、rank X 番目のステップで必要な直交基底が全て求まるためそこで処理を終える。
このとき残された列が、除外される説明変数に相当する（これらは、rank X 番目のステップまでに求めた基底で表現できてしまう）。
畳屋民也
summary で表示されるメッセージの表示条件
こちら で実行したように、入力説明変数に完全な多重共線性がある場合、summary 関数を実行した際に
Coefficients: (1 not defined because of singularities)
のような表示が現れる。
これは、下記で制御されている：
このうち、
df[3L] - df[1L]
は、（計画行列Xのカラム数） - （Xの階数）を表しており、Xが full-rank でない場合に表示処理が走るようになっている。
いか、それぞれどのように取得しているかを辿る。
df はどこから来ているか？
以下で定義されている：
df[1L] (p)
p は計画行列 X の rank と一致する
df[3L] (NCOL(Qr$qr))
Qr$qr は「QR分解後の上三角行列にXを再生成するために必要な情報を加えたもの」を表しているとのことだが、それがどうやら計画行列 X と同じ行・列数を持っているらしい。
てっきりXの階数（線形独立な成分の数）と一致するのかと思っていただけに、不思議。
参考： df[2L] (rdf)
残差の自由度 (
補足： Qr$qr の中身
以前実行した内容 における、 Qr$qr に相当する内容をみたところ、以下のようであった：
res$qr$qr
(Intercept) x1 x2
1 -2.2360680 -60.3738354 -1.811215e+02
2 0.4472136 23.0217289 6.906519e+01
3 0.4472136 0.5615152 1.820224e-14
4 0.4472136 -0.7416015 7.807201e-01
5 0.4472136 0.1271430 4.879500e-01
入力した x1, x2 の値は以下の通りなので、全く異なる値であることがわかる。
x1 <- c(30, 20, 15, 45, 25)
x2 <- x1 * 3
lm 関数で QR 分解時に行っている pivoting はなぜ入力データの列順に依存するか？
以前のメモ で lm 関数内で QR 分解を行う際に pivoting という列の入れ替えを行っていると述べた。
通常の pivoting ではノルムの大きさの順で列の並び替えを行うため、入力データの並び順を変えても出力結果は変わらない（補足の SciPy での例を参照）。
しかし、lm 関数では結果が入力値の列の並び順に依存する。
実際、説明変数を (x1, x2) の順で指定すると、x2 が除外されたのに対して、
下記のように順番を入れ替えて (x2, x1) にすると今度は x1 が除外される。
res2 <- lm(y~x2+x1, data=df)
summary(res2)
Call:
lm(formula = y ~ x2 + x1, data = df)
Residuals:
1 2 3 4 5
0.01590 -0.06889 0.01476 -0.01382 0.05205
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.0057122 0.0652887 30.72 7.58e-05 ***
x2 0.1017199 0.0007531 135.06 8.95e-07 ***
x1 NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05202 on 3 degrees of freedom
Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
F-statistic: 1.824e+04 on 1 and 3 DF, p-value: 8.949e-07
ではどのようなアルゴリズムで順番が決まるのかを見ていく。
lm 関数内で使われている QR 分解について
以前のメモ でも見たように、lm 関数では以下で QR 分解を行なっている：
[https://github.com/SurajGupta/r-source/blob/master/src/appl/dqrdc2.f]
これの冒頭を見ると、以下のように書かれている：
どうやら、pivoting を適用するカラムを限定しているようである。
同じく dqrdc2 を利用している qr 関数のドキュメントを見ると、以下のような記述があった：
qr() uses a modified version of LINPACK's DQRDC, called ‘dqrdc2’. It differs by using the tolerance tol for a pivoting strategy which moves columns with near-zero 2-norm to the right-hand edge of the x matrix. This strategy means that sequential one degree-of-freedom effects can be computed in a natural way.
https://stat.ethz.ch/R-manual/R-devel/library/base/html/qr.html
これによると、「ノルムが大きい順に並び替える」ではなく、先頭から見ていってもしノルムが0に近ければ列の最後に持っていく、のようにしていると考えられる。
なお、上記の dqrdc2 の実装は FORTRAN の線形代数ライブラリである LINPACK （NumPy や SciPy で使われている LAPACK の元となったもの）の QR 分解 routine dqrdc がベースとなっており、オリジナル版もこちらにある。
補足： SciPy の scipy.linalg.qr を使った場合
これがもし、SciPy のscipy.linalg.qr であれば、入力データの列の順番を入れ替えても x2→定数項→x1 の順に並ぶ。
import scipy as sp
x = np.column_stack([np.ones(df.shape[0]), df[["x1", "x2"]].to_numpy()])
sp.linalg.qr(x, pivoting=True, mode="economic")
(array([[-0.46429394, -0.03758018, 0.86117489],
[-0.30952929, -0.4434461 , -0.02171927],
[-0.23214697, -0.64637907, -0.16310449],
[-0.69644091, 0.57121871, -0.30227578],
[-0.38691162, -0.24051314, -0.37407535]]),
array([[-1.93842720e+02, -2.08932273e+00, -6.46142399e+01],
[ 0.00000000e+00, -7.96699780e-01, -3.36655700e-15],
[ 0.00000000e+00, 0.00000000e+00, 7.19549626e-15]]),
array([2, 0, 1], dtype=int32))
x_swap = np.column_stack([np.ones(df.shape[0]), df[["x2", "x1"]].to_numpy()])
sp.linalg.qr(x_swap, pivoting=True, mode="economic")
(array([[-0.46429394, -0.03758018, 0.86117489],
[-0.30952929, -0.4434461 , -0.02171927],
[-0.23214697, -0.64637907, -0.16310449],
[-0.69644091, 0.57121871, -0.30227578],
[-0.38691162, -0.24051314, -0.37407535]]),
array([[-1.93842720e+02, -2.08932273e+00, -6.46142399e+01],
[ 0.00000000e+00, -7.96699780e-01, -3.36655700e-15],
[ 0.00000000e+00, 0.00000000e+00, 7.19549626e-15]]),
array([1, 0, 2], dtype=int32))