在上次的代碼重寫中使用了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( ∑ (預測值 – 實際值)2 )，這就意味著如果人工計算的話就需要窮舉所有的函式，計算他們的損失然后找出最小的那個函式，
再來說說梯度下降，梯度下降是迭代法的一種，通過更新梯度來找到損失最小的那個函式，不知你發現問題了嗎？最小二乘法與梯度下降是在做同一件事情，也就是最優化問題，兩個是并行的關系，并不存在誰解決誰，百度百科中關于梯度下降的詞條中提到：“梯度下降是迭代法的一種，可以用于求解最小二乘問題(線性和非線性都可以)，在求解機器學習演算法的模型引數，即無約束優化問題時，梯度下降（Gradient Descent）是最常采用的方法之一，另一種常用的方法是最小二乘法，”這里很清晰的指出了最小二乘法和梯度下降法的關系，
再來說說另一個我之前并不知道的概念：最小二乘準則，百度百科中提到這是一種對于偏差程度的評估準則，與上兩者不同，上述的演算法都是基于最小二乘準則提出的對于最小二乘法優化問題的解決方案，也就是如果不窮舉的話如何找到最小二乘法的最優解，
列出我查閱的資料：
1、知乎-最小二乘法和梯度下降法有哪些區別？
2、百度百科-梯度下降
3、百度百科-最小二乘法
總結
1、雖然找到了sklearn.LinearRegression 類中對于線性回歸的演算法及實作，但發現并沒有使用到梯度下降法，而是使用最小二乘法找到最優解，解開了我對最小二乘法與梯度下降到誤解，但由于之前并未從事過演算法研究與數學分析，對相應的演算法一知半解，所以這里的代碼難以看懂，只能就此作罷，學習了相應的演算法之后再來學習代碼實作， 2、在學習原始碼的程序中以及寫這篇文章的程序中發現對于python的有些概念還是不太清晰，比如函式和方法還有fun(x)(y) 的呼叫方式，所以能看到上文中有些使用“函式”有些使用“方法”，可能并不對應，但之后熟悉了才能修改， 3、原始碼中的注釋充斥著許多數學詞匯，讀起來讓我例外頭疼，幾乎都要使用翻譯軟體才能理解，同時有些平常使用的詞匯我也不懂，這個時候英語的作用就十分必要了， 4、總的來說，基礎不扎實，水平不高，所以對于其中的精髓難以理解，同時可能文章中錯漏百出，但我并未發現，這就是目前的問題或者說困境，勤加學習才能脫離，轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/281150.html
標籤：其他
下一篇：frida rpc遠程呼叫