深度學習03-sklearn.LinearRegression 原始碼學習

RabbitKeeper發表於2021-04-28
深度學習原始碼

在上次的程式碼重寫中使用了sklearn.LinearRegression 類進行了線性迴歸之後猜測其使用的是常用的梯度下降+反向傳播演算法實現，所以今天來學習它的原始碼實現。但是在看到原始碼的一瞬間突然有種懷疑人生的感覺，我是誰？我在哪？果然大佬的程式碼只能讓我膜拜。

在一目十行地看完程式碼之後，我發現了一個問題，梯度的單詞是gradient，一般在程式碼中會使用縮寫grad 來表示梯度，而在這個程式碼中除了Gram 之外竟然沒有一個以'g' 開頭的單詞，更不用說gradient 了。那麼程式碼中包括註釋壓根沒提到過梯度，是不是說明這裡根本沒有使用梯度下降演算法呢，換言之就是是否還有其他方法來實現最小二乘法的線性迴歸呢？帶著這個疑問，我開始仔細閱讀LinearRegression.fit() 函式。

首先都是引數處理，這裡雖然看不太懂但是也能大概知道他在做什麼，所以可以跳過。

然後來到了核心程式碼，核心程式碼中使用的幾個判斷：

 1 if self.positive:
 2     if y.ndim < 2:
 3         pass
 4     else:
 5         pass
 6 elif sp.issparse(X):
 7     if y.ndim < 2:
 8         pass
 9     else:
10         pass
11 else:
12     pass

self.positive 是在使用密集矩陣的時候設定的引數，y.ndim 表示y 的維度，簡單來說就是y 中有幾個[]，所以大概能知道程式碼將密集矩陣與稀疏矩陣區分開，並且將一維矩陣與多維矩陣區分開，意味著不同的類別使用不同的方法。

if self.positive 分段解析

1 if self.positive:
2     if y.ndim < 2:
3         self.coef_, self._residues = optimize.nnls(X, y)
4     else:
5         # scipy.optimize.nnls cannot handle y with shape (M, K)
6         outs = Parallel(n_jobs=n_jobs_)(
7             delayed(optimize.nnls)(X, y[:, j])
8             for j in range(y.shape[1]))
9         self.coef_, self._residues = map(np.vstack, zip(*outs))

可以看出y 的維度小於2 的話使用optimize.nnls() 方法，否則進行其他處理，因為“scipy.optimize.nnls cannot handle y with shape (M, K)”，但看到之後也呼叫了optimize.nnls，所以應該是將矩陣處理成可以使用的樣子。

並且值得注意的是這裡使用了Parallel(n_jobs=n_jobs_)(delayed(optimize.nnls)(X, y[:, j])for j in range(y.shape[1])) 的呼叫方式，也就是形如fun(x)(y) 的方式，這意味著函式內定義了另一個函式，第一個括號是fun 的引數，第二個括號是給fun 函式內定義的函式的引數。

elif sp.issparse(X) 分段解析

 1 elif sp.issparse(X):
 2 X_offset_scale = X_offset / X_scale
 3 
 4 
 5 def matvec(b):
 6     return X.dot(b) - b.dot(X_offset_scale)
 7 
 8 
 9 def rmatvec(b):
10     return X.T.dot(b) - X_offset_scale * np.sum(b)
11 
12 
13 X_centered = sparse.linalg.LinearOperator(shape=X.shape,
14                                           matvec=matvec,
15                                           rmatvec=rmatvec)
16 
17 if y.ndim < 2:
18     out = sparse_lsqr(X_centered, y)
19     self.coef_ = out[0]
20     self._residues = out[3]
21 else:
22     # sparse_lstsq cannot handle y with shape (M, K)
23     outs = Parallel(n_jobs=n_jobs_)(
24         delayed(sparse_lsqr)(X_centered, y[:, j].ravel())
25         for j in range(y.shape[1]))
26     self.coef_ = np.vstack([out[0] for out in outs])
27     self._residues = np.vstack([out[3] for out in outs])

可以看到先是對資料進行了處理，然後呼叫了sparse_lsqr() 函式。

剩餘分段解析

1 else:
2 self.coef_, self._residues, self.rank_, self.singular_ = \
3     linalg.lstsq(X, y)
4 self.coef_ = self.coef_.T
5 
6 if y.ndim == 1:
7     self.coef_ = np.ravel(self.coef_)
8 self._set_intercept(X_offset, y_offset, X_scale)

使用了linalg.lstsq() 函式。

以上我們可以看到程式碼中一共使用了3 個方法來實現線性迴歸：optimize.nnls()、sparse_lsqr()、linalg.lstsq()

optimize.nnls() 分析

NNLS 即非負正則化最小二乘法，程式碼實現由scipy.optimize.nnls 提供，這裡只是將其封裝起來，在原始碼的註釋中提到該演算法的FORTRAN 程式碼在Charles L. Lawson 與Richard J. Hanson 兩位教授於1987 年所著的《Solving Least Squares Problems》中釋出。“The algorithm is an active set method. It solves the KKT (Karush-Kuhn-Tucker) conditions for the non-negative least squares problem.” 可惜由於本人水平有限，並不能從書中或者此處的程式碼中學到該演算法的精髓，只能先挖一個坑，以後有所提高了再來研究該演算法。

sparse_lsqr() 分析

LSQR 即最小二乘QR分解演算法，程式碼實現由scipy.sparse.linalg.lsqr 類提供，這裡只是將其封裝起來，在文件中可以看到：

LSQR uses an iterative method to approximate the solution.  The number of iterations required to reach a certain accuracy depends strongly on the scaling of the problem.  Poor scaling of the rows or columns of A should therefore be avoided where possible.

 同時也給出了參考文獻：

[1] C. C. Paige and M. A. Saunders (1982a). "LSQR: An algorithm for sparse linear equations and sparse least squares", ACM TOMS 8(1), 43-71.
[2] C. C. Paige and M. A. Saunders (1982b). "Algorithm 583.  LSQR: Sparse linear equations and least squares problems", ACM TOMS 8(2), 195-209.
[3] M. A. Saunders (1995).  "Solution of sparse rectangular systems using LSQR and CRAIG", BIT 35, 588-604.

 可以看到LSQR 演算法是Paige 和Saunders 於1982 年提出的一種方法，但水平有限，暫時並不清楚其中原理。

linalg.lstsq() 分析

LSTSQ 是 LeaST SQuare （最小二乘）的意思，也就是普通的最小二乘法。程式碼實現由scipy.linalg.lstsq 提供，這裡只是將其封裝起來。通過官方文件提供的原始碼連結，我找到了lstsq 函式的原始碼，註釋中提到了：

Which LAPACK driver is used to solve the least-squares problem.
Options are ``'gelsd'``, ``'gelsy'``, ``'gelss'``. Default(``'gelsd'``) is a good choice.  However, ``'gelsy'`` can be slightly faster on many problems.  ``'gelss'`` was used historically.  It is generally slow but uses less memory.

 也就是說有三個選項：gelsd(預設推薦)、gelsy(可能稍快)、gelss(使用記憶體少)，那麼來看看他們分別使用什麼方法來解決最小二乘法吧。

 1 if driver in ('gelss', 'gelsd'):
 2     if driver == 'gelss':
 3         lwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond)
 4         v, x, s, rank, work, info = lapack_func(a1, b1, cond, lwork,
 5                                                 overwrite_a=overwrite_a,
 6                                                 overwrite_b=overwrite_b)
 7 
 8     elif driver == 'gelsd':
 9         if real_data:
10             lwork, iwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond)
11             x, s, rank, info = lapack_func(a1, b1, lwork,
12                                            iwork, cond, False, False)
13         else:  # complex data
14             lwork, rwork, iwork = _compute_lwork(lapack_lwork, m, n,
15                                                  nrhs, cond)
16             x, s, rank, info = lapack_func(a1, b1, lwork, rwork, iwork,
17                                            cond, False, False)
18 elif driver == 'gelsy':
19         lwork = _compute_lwork(lapack_lwork, m, n, nrhs, cond)
20         jptv = np.zeros((a1.shape[1], 1), dtype=np.int32)
21         v, x, j, rank, info = lapack_func(a1, b1, jptv, cond,
22                                           lwork, False, False)

 

 看到呼叫了lapack_func，但是找了一下spicy 並沒有發現這個函式，於是搜尋lapack_func 找到了定義：

lapack_func, lapack_lwork = get_lapack_funcs((driver,'%s_lwork' % driver),(a1, b1))

原來它呼叫了get_lapack_funcs 函式，於是去找該函式的文件，從文件中看到了該方法的用處：“This routine automatically chooses between Fortran/C interfaces. Fortran code is used whenever possible for arrays with column major order. In all other cases, C code is preferred.”，原來這是一個Fortran 語言和C 語言的介面且首選C 語言，並且從原始碼中我看到了它是呼叫_get_funcs 實現，那3 個單詞是C 語言的函式名，所幸我對C 語言的熟悉程度打過Python，於是去找C 語言API 進行學習。

gelsd：Computes the minimum-norm solution to a linear least squares problem using the singular value decomposition of A and a divide and conquer method. 分治法的奇異值分解找出最優解。

gelss：Computes the minimum-norm solution to a linear least squares problem using the singular value decomposition of A. 使用奇異值分解找出最優解。

gelsy：Computes the minimum-norm solution to a linear least squares problem using a complete orthogonal factorization of A. 使用完全正交分解的方式找出最優解。

至此對於所有函式的分析基本上是結束了，密集矩陣使用的是NNLS 演算法，稀疏矩陣使用的是LSQR 演算法，其他的使用的則是最常用的最小二乘法演算法。回到開頭的問題，為什麼沒有使用梯度下降呢？難道梯度下降並不是解決最小二乘或者說是線性迴歸的演算法嗎？在網上查閱了很多材料之後我發現自己的問題：我陷入了誤區。

先來說說最小二乘法，最小二乘法在我初高中的時候學習簡單線性迴歸的時候就接觸了，根據百度百科的詞條，最小二乘法（又稱最小平方法）是一種數學優化技術。它通過最小化誤差的平方和尋找資料的最佳函式匹配。利用最小二乘法可以簡便地求得未知的資料，並使得這些求得的資料與實際資料之間誤差的平方和為最小，也就是說最小二乘法的公式是目標函式 = MIN( ∑ (預測值 – 實際值)² )。這就意味著如果人工計算的話就需要窮舉所有的函式，計算他們的損失然後找出最小的那個函式。

再來說說梯度下降，梯度下降是迭代法的一種，通過更新梯度來找到損失最小的那個函式，不知你發現問題了嗎？最小二乘法與梯度下降是在做同一件事情，也就是最優化問題，兩個是並行的關係，並不存在誰解決誰。百度百科中關於梯度下降的詞條中提到：“梯度下降是迭代法的一種，可以用於求解最小二乘問題(線性和非線性都可以)。在求解機器學習演算法的模型引數，即無約束優化問題時，梯度下降（Gradient Descent）是最常採用的方法之一，另一種常用的方法是最小二乘法。”這裡很清晰的指出了最小二乘法和梯度下降法的關係。

再來說說另一個我之前並不知道的概念：最小二乘準則。百度百科中提到這是一種對於偏差程度的評估準則，與上兩者不同，上述的演算法都是基於最小二乘準則提出的對於最小二乘法優化問題的解決方案，也就是如果不窮舉的話如何找到最小二乘法的最優解。

列出我查閱的資料：

1、知乎-最小二乘法和梯度下降法有哪些區別？

2、百度百科-梯度下降

3、百度百科-最小二乘法

總結

1、雖然找到了sklearn.LinearRegression 類中對於線性迴歸的演算法及實現，但發現並沒有使用到梯度下降法，而是使用最小二乘法找到最優解，解開了我對最小二乘法與梯度下降到誤解，但由於之前並未從事過演算法研究與數學分析，對相應的演算法一知半解，所以這裡的程式碼難以看懂，只能就此作罷，學習了相應的演算法之後再來學習程式碼實現。
2、在學習原始碼的過程中以及寫這篇文章的過程中發現對於python的有些概念還是不太清晰，比如函式和方法還有fun(x)(y) 的呼叫方式，所以能看到上文中有些使用“函式”有些使用“方法”，可能並不對應，但之後熟悉了才能修改。
3、原始碼中的註釋充斥著許多數學詞彙，讀起來讓我異常頭疼，幾乎都要使用翻譯軟體才能理解，同時有些平常使用的詞彙我也不懂，這個時候英語的作用就十分必要了。
4、總的來說，基礎不紮實，水平不高，所以對於其中的精髓難以理解，同時可能文章中錯漏百出，但我並未發現，這就是目前的問題或者說困境，勤加學習才能脫離。

相關文章