R の lm 関数

結果

default では、注意書き付きで変数選択したうえで結果を出力した。
singular.ok = FALSE をオプション指定すると、エラーが出る。

実行コード

x2 = x1 * 3 という多重共線性のあるデータを入力。

x1 <- c(30, 20, 15, 45, 25)
x2 <- x1 * 3
y <- 2 + 0.3 * x1 + 0.1 * rnorm(5)

df <- data.frame(x1, x2, y)

res <- lm(y~x1+x2, data=df)
summary(res)

出力結果

x2 を除外した推定結果を返している。

Call:
lm(formula = y ~ x1 + x2, data = df)

Residuals:
       1        2        3        4        5 
 0.00140 -0.09743  0.07571  0.01315  0.00718 

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.062194   0.090078   22.89 0.000183 ***
x1          0.296699   0.003117   95.18 2.56e-06 ***
x2                NA         NA      NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07177 on 3 degrees of freedom
Multiple R-squared:  0.9997,	Adjusted R-squared:  0.9996 
F-statistic:  9059 on 1 and 3 DF,  p-value: 2.557e-06

`singular.ok = FALSE` を指定した場合

実行コード

res <- lm(y~x1+x2, data=df, singular.ok = FALSE)

出力結果

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
  singular fit encountered

畳屋民也

statsmodels の OLS (Python)

いずれの場合も summary に注意書き付きで実行結果が返ってくる。
default のオプションでは Moore-Penrose の逆行列を用いた推定値が返ってくる
method="qr" を指定することで QR 分解を用いた推定値が返ってくる。
- R の lm 関数の時と異なり、変数選択は行わない。
- 推定が不安定なのか、係数の大きさが極端な推定値が返ってきたことに加え、係数標準誤差がうまく計算できなかった。

実行結果

x2 = x1 * 3 という完全な多重共線性のあるデータを入力し、最小二乗法によりモデル推定を行った。

モデル作成

default のオプションでは、Moore-Penrose の逆行列(pinv)を用いた推定を行う。

import pandas as pd
import numpy as np
import statsmodels.api as sm

np.random.seed(0)

# データ準備
df = pd.DataFrame({"x1": [30, 20, 15, 45, 25]})
df["x2"] = df["x1"] * 3
df["y"] = 2 + 0.3 * df["x1"] +  0.1 * np.random.randn(df.shape[0])

# モデル作成（推定方法： Moore-Penrose の逆行列）
mod = sm.OLS(df["y"], sm.add_constant(df[["x1", "x2"]]))
res = mod.fit()

推定結果

特に問題なく推定結果がおこなわれる。
固有値が小さいため summary に注意書きは表示される。

print(res.summary())

        OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 1.824e+04
Date:                Tue, 13 Feb 2024   Prob (F-statistic):           8.95e-07
Time:                        14:27:43   Log-Likelihood:                 8.9634
No. Observations:                   5   AIC:                            -13.93
Df Residuals:                       3   BIC:                            -14.71
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.0057      0.065     30.721      0.000       1.798       2.213
x1             0.0305      0.000    135.062      0.000       0.030       0.031
x2             0.0915      0.001    135.062      0.000       0.089       0.094
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   2.383
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.365
Skew:                          -0.555   Prob(JB):                        0.833
Kurtosis:                       2.278   Cond. No.                     1.44e+16
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 2.02e-28. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
/usr/local/lib/python3.10/dist-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.
  warn("omni_normtest is not valid with less than 8 observations; %i "

QR分解を用いた最小二乗法による推定結果

fit 時のオプションで、R の lm 関数と同様の QR 分解を用いた推定を行うよう指定してみた。

pinv と同様に summary に注意書き付きで実行結果が返ってくるものの、
係数推定値が $-2.622 \times 10^{12}, \, 8.741 \times 10^{11}$ と極端な値となっていることに加えて、
係数分散が負として計算されたため標準誤差が nan になった。

res_qr = mod.fit(method="qr")
print(res_qr.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       1.000
Model:                            OLS   Adj. R-squared:                  1.000
Method:                 Least Squares   F-statistic:                 1.400e+04
Date:                Tue, 13 Feb 2024   Prob (F-statistic):           1.33e-06
Time:                        14:26:42   Log-Likelihood:                 8.3018
No. Observations:                   5   AIC:                            -12.60
Df Residuals:                       3   BIC:                            -13.38
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.9985      0.075     26.816      0.000       1.761       2.236
x1         -2.622e+12        nan        nan        nan         nan         nan
x2          8.741e+11        nan        nan        nan         nan         nan
==============================================================================
Omnibus:                          nan   Durbin-Watson:                   1.447
Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.287
Skew:                          -0.536   Prob(JB):                        0.866
Kurtosis:                       2.522   Cond. No.                     1.44e+16
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 2.02e-28. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.
/usr/local/lib/python3.10/dist-packages/statsmodels/stats/stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.
  warn("omni_normtest is not valid with less than 8 observations; %i "
/usr/local/lib/python3.10/dist-packages/statsmodels/regression/linear_model.py:1884: RuntimeWarning: invalid value encountered in sqrt
  return np.sqrt(np.diag(self.cov_params()))

畳屋民也

scikit-learn の LinearRegression (Python)

特に Warning や注意書き等の表示なく結果が返ってくる
推定された切片・回帰係数の値は、statsmodels で Moore-Penrose 逆行列 (pinv) により求めたものと同じ値

データ準備

前回と同様に説明変数同士に x2 = 3*x1 という関係を持たせ、目的変数は y = 2 + 0.3 * x1 + eps のように生成した。

import pandas as pd
import numpy as np

np.random.seed(0)

# データ準備
df = pd.DataFrame({"x1": [30, 20, 15, 45, 25]})
df["x2"] = df["x1"] * 3
df["y"] = 2 + 0.3 * df["x1"] +  0.1 * np.random.randn(df.shape[0])

モデル作成

以下のように Warning や注意書き等が出ることなく結果が返ってきた。

from sklearn.linear_model import LinearRegression

reg = LinearRegression()
reg.fit(df[["x1", "x2"]], df["y"])

print(reg.intercept_)
print(reg.coef_)

2.005712241700081
[0.03051598 0.09154795]

これらの推定値は、statsmodels によるモデルで default オプション method="pinv" を使った際の結果と一致している。

畳屋民也

scikit-learn の LinearRegression の内部の最適化アルゴリズムについて

scikit-learn の LinearRegression で最小二乗法によりモデル推定を行う際に、どのような最適化アルゴリズムを用いているかを調べた。

上記ドキュメントにあるように、default では SciPy の scipy.linalg.lstsq を使っていた。
（positive=True として係数推定値が正の値を取るように制約を加えた場合、scipy.optimize.nnls を使う。）

From the implementation point of view, this is just plain Ordinary Least Squares (scipy.linalg.lstsq) or Non Negative Least Squares (scipy.optimize.nnls) wrapped as a predictor object.

さらに辿ると、LAPACK の関数（を wrap したもの）を呼び出しており、Moore-Penrose の逆行列を用いたのと同等のノルム最小推定量を求めていた（ように読めた）。

scikit-learn のコードの該当箇所

デフォルトオプションの positive=False かつ入力が疎行列でなければ、以下の箇所で scipy.linalg.lstsq が呼び出される。

scipy.linalg.lstsq の中身

default では、入力値が実数なら scipy.linalg.lapack.dgelsd が呼び出される。
これは、LAPACK (Linear Algebra PACKage) の dgelsd 関数を wrap したものだという。

dgelsd とは？

DGELSD computes the minimum-norm solution to a real linear least
squares problem:
minimize 2-norm(| b - A*x |)
using the singular value decomposition (SVD) of A. A is an M-by-N
matrix which may be rank-deficient.

これによると、特異値分解を用いて最小二乗法の最小ノルム解を計算するという。

前掲のドキュメントに基づくと、より詳細には、

計画行列 X を Householder 変換により bidiagonal form （二重対角行列？）にして、"bidiagonal least squares problem (BLS)"（二重対角最小二乗問題？）に落とし込む
分割統治法により BLS を解く
Householder 変換の逆変換をかける

ということを行っているらしい。

bidiagonal form について

dgelsd を辿る

実質的に、Moore-Penrose の逆行列（擬逆行列）を用いた最小ノルム推定量を求めていると見て良さそうである。

ソースコードを見ると条件訳して処理が分岐しているものの、最終的には計画行列 X を bidiagonal form に変換したのち、dlalsd という別の関数を適用している。

こちらをみると、以下のような操作で $||Ax - B||_2$ を最小化する $x$ を求めているように読み取れる（この場合、 $A$ が bidiagonal form にした計画行列、 $x$ が回帰係数、 $B$ が目的変数に相当）：

dlasda という関数で bidiagonal 行列 $A$ を特異値分解する
- $A = U \Sigma V^\top$ （ $\Sigma$ は対角行列で、対角成分は特異値）

行列 $B$ と右特異行列 $U$ をかける
- $U^\top B$

特異値の逆数をかける（ただし、特異値が一定値以下の場合は0とみなし除外）
- $\Sigma^{-1} U^\top B$

左から右特異行列をかける
- $V \Sigma^{-1} U^\top B$

なお、特異値分解と Moore-Penrose 逆行列の関係は以下を参照：

参考： bidiagonal least squares problem (BLS) について

畳屋民也

statsmodels で QR 分解を選んだ時の処理を辿る

こちらで見たように、statsmodels の OLS で fit の際に method="qr" を指定すると最適化アルゴリズムとして QR 分解が使用される。

しかし、同じQR分解を用いる R の lm 関数を使用した場合とは異なる結果が得られた。

この理由を調べるため、statsmodels の線形回帰のコードを辿った。

statsmodels で QR 分解が行われている箇所

statsmodels の OLS クラスはこちらで定義されている：

さらに元を辿ると、fit メソッドは以下の RegressionModel から継承している：

method="qr" を選択すると、以下の処理が行われる：

この中のL348で、NumPy の numpy.linalg.qr 関数を呼び出している。

NumPy による QR 分解

NumPy の numpy.linalg.qr 関数の中では LAPACK の QR 分解の routine を呼び出しているらしい。

`numpy.linalg.qr` を使ってみる

計画行列 $X$ を実際に numpy.linalg.qr 関数で QR 分解してみた。

x = np.column_stack([np.ones(df.shape[0]), df[["x1", "x2"]].to_numpy()])
Q, R = np.linalg.qr(x, mode="reduced")
Q, R

(array([[-0.4472136 ,  0.13031167, -0.85065652],
        [-0.4472136 , -0.30406057,  0.40047151],
        [-0.4472136 , -0.52124669,  0.01420386],
        [-0.4472136 ,  0.78187004,  0.31988394],
        [-0.4472136 , -0.08687445,  0.11609722]]),
 array([[-2.23606798e+00, -6.03738354e+01, -1.81121506e+02],
        [ 0.00000000e+00,  2.30217289e+01,  6.90651866e+01],
        [ 0.00000000e+00,  0.00000000e+00, -4.49386684e-14]]))

この場合、直交行列 $Q$ は $5 \times 3$ 、上三角行列 $R$ は $3\times 3$ で返ってきている。
行列 R の rank は計画行列 $X$ と同様に2であり、正則でない。

np.linalg.matrix_rank(R)

R言語のlm関数であれば、変数選択により $Q \in \mathcal{R}^{5\times 2}, \, R \in \mathcal{R}^{2 \times 2}$ 相当の結果を返していたところである。

次元を減らして回帰係数を計算してみる

ここで、R言語のlm関数を参考に直交行列 $Q$ は2列目まで、上三角行列 $R$ は先頭の $2 \times 2$ 主部分行列だけ取り出して回帰係数を算出して見た。
（ただし後述するが、lm 関数の場合と異なり numpy.linalg.qr では pivoting は行っていない。）

effects2 = np.dot(Q[:, :2].T, df["y"])
beta2 = np.linalg.solve(R[:2, :2], effects2)
beta2

array([2.00571224, 0.30515984])

この結果は、Rのlm関数を使った場合とほぼ一致する。
以前の実行結果と僅かな差はあるものの、目的変数 $y$ を生成する際の乱数の違いによるものである。
同一のデータを用いてあらためて lm 関数でモデル推定を行った結果、同一の係数推定値が得られた。

`numpy.linalg.qr` では pivoting を行っているか？

QRアルゴリズムでは、pivoting といって基底を計算する iteration の際になるべくそれまで計算した基底と直行した列を選ぶように列の並び替えを行うことがある。

R の lm では pivoting を（部分的に）行っていたが、 numpy.linalg.qr ではどうかを調べてみた。

結論から言うと、numpy.linalg.qr では pivoting は行っていない。

numpy.linalg.qr では、以下の箇所で dgeqrf （複素数の場合は zgdeqrf）という LAPACK の関数を呼び出している。^[1]

dgeqrf のドキュメントには以下のような記述がある：

F08AEF (DGEQRF) forms the QR factorization of an arbitrary rectangular real m by n matrix. No
pivoting is performed.
https://support.nag.com/numeric/fl/nagdoc_fl26/pdf/f08/f08aef.pdf

ということで、pivoting は行っていないらしい。

なお、 SciPy の QR 分解関数 scipy.linalg.qr では pivoting を行うオプションが存在する。

脚注

正確には、その後でもう一箇所で orgqrf という LAPACK の関数を呼び出している。筆者の理解では、geqrf で返ってきた値（前述のドキュメントでは The matrix Q is not formed explicitly but is represented as a product of min(m, n) elementary reflectors とのこと）をもとに、orgqrf で直交行列 Q を構成している。 ↩︎

畳屋民也

補足： SciPy の QR 分解（pivoting あり）

SciPy の scipy.linalg.qr では、pivoting=True を指定すると pivoting ありの QR 分解を行う。

これがどこで行われているかを辿った。

`scipy.linalg.qr` のコード

以下で scipy.linalg.qr が定義されている。

[https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/linalg/_decomp_qr.py#L26-L176]

このなかで、pivoting=True の場合は以下のように get_lapack_funcs 関数によって LAPACK の geqp3 という関数が呼び出される。

LAPACK の `dgeqp3`

入力が実数行列の場合は、以下の dgeqp3 が呼び出される：

[https://github.com/Reference-LAPACK/lapack/blob/c6bc40164ace127bc308da37b135791a9ecb4667/SRC/dgeqp3.f]

この中で実際に pivoting しながら行列分解を行なっている箇所は、以下で呼び出している DLAQPS という routine と見られる：

[https://github.com/Reference-LAPACK/lapack/blob/c6bc40164ace127bc308da37b135791a9ecb4667/SRC/dlaqps.f]

R の lm 関数

実行コード

出力結果

singular.ok = FALSE を指定した場合

statsmodels の OLS (Python)

実行結果

モデル作成

推定結果

QR分解を用いた最小二乗法による推定結果

scikit-learn の LinearRegression (Python)

データ準備

モデル作成

scikit-learn の LinearRegression の内部の最適化アルゴリズムについて

scikit-learn のコードの該当箇所

scipy.linalg.lstsq の中身

dgelsd とは？

dgelsd を辿る

参考： bidiagonal least squares problem (BLS) について

statsmodels で QR 分解を選んだ時の処理を辿る

statsmodels で QR 分解が行われている箇所

NumPy による QR 分解

numpy.linalg.qr を使ってみる

次元を減らして回帰係数を計算してみる

numpy.linalg.qr では pivoting を行っているか？

補足： SciPy の QR 分解（pivoting あり）

scipy.linalg.qr のコード

LAPACK の dgeqp3

`singular.ok = FALSE` を指定した場合

`numpy.linalg.qr` を使ってみる

`numpy.linalg.qr` では pivoting を行っているか？

`scipy.linalg.qr` のコード

LAPACK の `dgeqp3`