常用算法分析——最小二乘法

1甥雕、引言

最小二乘法應該是我們最早接觸的一種數值估計算法。它的特殊形式——一元線性回歸服协，被廣泛地應用于多種數值統(tǒng)計分析場合。例如，在驗證歐姆定律（U = IR）時，通常的實驗方法是分別測量出多個不同電壓U_i下恰响，通過電阻的電流值I_i，然后將這些(U_i, I_i)觀測點涌献，代入到一元最小二乘公式(1-1)中胚宦，便可計算出 $\hat{R}$ 。

$\begin{cases}a&=&\frac{\sum{xy}-\frac{1}{N}\sum{x}\sum{y}}{\sum{x^2}-\frac{1}{N}(\sum{x})^2}\\ b&=&\frac{1}{N}\sum{y}-\frac{a}{N}\sum{x}\end{cases}$ (1-1)

由此可得出線性擬合式(1-2)燕垃，

$\hat{y}=a\hat{x}+b$ (1-2)

其中枢劝， $\hat{y}=\hat{U},\ \hat{x}=\hat{I},\ a=\hat{R},\ b$ 是殘差。通過此方法將觀測點及擬合曲線繪制在同一個直角坐標系中卜壕，正常情況下可以直觀地看到您旁，觀測點會均勻分布在直線附近，且每個點的殘差平方和（即方差）最小轴捎。
“最小二乘法”由此得名鹤盒。

2、普通最小二乘法（OLS）

最小二乘法顯然不只是一元線性回歸那么簡單侦副，它還可以應用于多元參數的擬合昨悼。本節(jié)將對普通最小二乘法（Ordinary Least Squares）的原理進行簡單的推導和證明。

2.1跃洛、高斯—馬爾可夫定理

高斯—馬爾可夫定理（the Gauss–Markov theorem，簡稱G-M定理）　在給定經典線性回歸的假定下终议，最小二乘估計量是具有最小方差的線性無偏估計量（即Best Linear Unbiased Estimator汇竭，簡稱BLUE)。

G-M定理共對OLS普通線性方程提出5個假設：
假設1（線性關系）：要求所有的母集團參數（population parameters）為常數穴张，用來保證模型為線性關系卿吐。即脂凶，如果母集團方程為 $a,b_1,b_2,...,b_k$ 必須為常數，ε為無法檢測的誤差項，即實驗過程中模型沒有包含的因素物赶。
假設2（隨機采樣）：假設我們有n個調查的樣本，那么這n個樣本必須是從母集團里面隨機抽樣得出的苦囱。
假設3（無完全共線關系）：在樣本（母集團）中栈暇，沒有獨立變量是常數，且獨立變量之間不能有完全共線性渐夸。
假設4（零條件均值）：母集團方程的誤差項均值為0嗤锉，且均值不受到獨立變量的影響，即E(ε | X₁, X₂, …, X_k)＝0墓塌。
假設5（同方差性）：誤差項ε的方差為一個固定不變的常數瘟忱，且不受到獨立變量的影響奥额，即Var(ε | X₁, X₂, …, X_k)＝σ²或Var(ε)＝σ²I。

G-M定理的意義在于访诱，當上述假定成立時垫挨，我們不需要再去尋找其它無偏估計量，因為沒有一個會優(yōu)于普通最小二乘估計量触菜。換言之九榔，如果存在一個好的線性無偏估計量，那么這個估計量的方差最多與普通最小二乘估計量的方差一樣小玫氢，而不會小于普通最小二乘估計量的方差帚屉。

2.2、最小二乘估計量推導

假設使用多輸入單輸出（MISO）線性系統(tǒng)對未知系統(tǒng)進行估計漾峡，則估計系統(tǒng)的輸入輸出應滿足

$y=\hat{\beta}_0+\hat{\beta}_1x_1+\hat{\beta}_2x_2+...+\hat{\beta}_kx_k=\sum_{j=0}^k\hat{\beta}_jx_j\ \ \ (x_0=1)$ (2-1)

此式稱為系統(tǒng)的估計量攻旦。接下來，定義第 i 次觀測值 y_i 應滿足

$y_i=\beta_0+\beta_1x_{1i}+\beta_2x_{2i}+...+\beta_kx_{ki}+\varepsilon_i$ (2-2)

其中生逸， $\small\beta_k$ 的估計量牢屋，ε_i 是第 i 次觀測時無法預測的誤差項。
由此可得

$\begin{array}{rcl}D&=& \sum_{i=1}^n \varepsilon_i^2\\ &=&\sum_{i=1}^n(y_i-\hat{y}_i)^2\\&=&\sum_{i=1}^n(y_i-\sum_{j=0}^k\hat{\beta}_jx_{ji})^2\end{array}$ (2-3)

若想求得使得式(2-3)取得極小值槽袄，那么需要對D分別求出 $\small\hat{\beta}_r$ 的偏導數烙无，并令其等于零，即遍尺，

$\frac{\partial D}{\partial\beta_r} =2\sum_{i=1}^{n}(y_i-\sum_{j=1}^k\hat{\beta}_jx_{ji})x_{ri}=2\sum_{i=1}^{n}x_{ri}y_i-2\sum_{i=1}^{n}(x_{ri}\sum_{j=0}^k\hat{\beta}_jx_{ji})=0$
$\Rightarrow \sum_{i=1}^{n}x_{ri}y_i=\sum_{i=1}^{n}(x_{ri}\sum_{j=0}^k\hat{\beta}_jx_{ji})$ (2-4)

其中截酷，r＝1, 2, 3, …, k，聯(lián)立這k個關于 $\hat{\beta}_r$ 的方程乾戏，即為法方程組（或正規(guī)方程組）迂苛，而法方程組的解即為式(2-1)的系數，通過此方法得到的估計量則稱為系統(tǒng)的最小二乘估計量鼓择∪茫可以證明，對于滿足上述5個假設的OLS方程組呐能，存在唯一實數解念搬。
特殊地，當k＝1時摆出，MISO退化成單輸入單輸出（SISO）的一元系統(tǒng)朗徊，則其法方程組可寫作：

$\left\{ \begin{array}{rcl} \sum_{i=1}^n(\hat{\beta}_0+\hat{\beta}_1x_i)&=& \sum_{i=1}^ny_i \\ \sum_{i=1}^n(\hat{\beta}_0x_i+\hat{\beta}_1x_i^2)&=&\sum_{i=1}^nx_iy_i \end{array}\right. \\ \Rightarrow \left\{ \begin{array}{rcl} n\hat\beta_0+(\sum_{i=1}^nx_i)\hat{\beta}_1&=& \sum_{i=1}^ny_i \\ (\sum_{i=1}^nx_i)\hat{\beta}_0+(\sum_{i=1}^nx_i^2)\hat{\beta}_1&=&\sum_{i=1}^nx_iy_i \end{array}\right.$

求解后，可得出式(1-1)的結果偎漫。

2.3荣倾、矩陣形式

上節(jié)已推導出計算最小二乘估計量的一般方法，但顯然骑丸，式(2-4)的形式過于繁瑣舌仍。對于多元線性回歸的場合妒貌，法方程組將會變得異常復雜，難以求解铸豁。同時灌曙，這種代數形式的公式也不便于得出更為一般的規(guī)律。因此节芥，本節(jié)將用矩陣的形式在刺，再次推導OLS的計算過程，并且嘗試討論法方程組解的存在問題头镊。
首先蚣驼，定義估計量

$\hat{y}={x}^\text{T}\hat{\boldsymbol{\beta}}$ (2-5)

其中， $x=(1,x_1,x_2,...,x_k)^\text{T},\ \hat{\boldsymbol{\beta}}=(\hat{\beta}_0,\hat{\beta}_1,\hat{\beta}_2,...,\hat{\beta}_k)^\text{T}$ 相艇。
誤差

$\boldsymbol{\varepsilon}=\boldsymbol{Y}-\hat{\boldsymbol{Y}}=\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}$ (2-6)

其中颖杏，
$\boldsymbol{\varepsilon}=(\varepsilon_1,\varepsilon_2,...,\varepsilon_n)^\text{T}$ ， $\boldsymbol{Y}=(y_1,y_2,...,y_n)^\text{T}$ 坛芽， $\boldsymbol{X}=(\boldsymbol{x}_1,\boldsymbol{x}_2,...,\boldsymbol{x}_n)^\text{T}=\left[\begin{matrix} 1&x_{11} &\cdots& x_{k1}\\1&x_{12}&&x_{k2}\\ \vdots&&\ddots&\vdots\\ 1&x_{1n} &\cdots& x_{kn}\end{matrix}\right]$ 留储。
取誤差向量的2-范數^[1]的平方

$\begin{array}{rcl} \boldsymbol{D}&=&\left\|\boldsymbol{\varepsilon}\right\|_2^2=\boldsymbol{\varepsilon}^\text{T}\boldsymbol{\varepsilon} \\ &=& (\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}})^\text{T}(\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}) \\ &=& (\boldsymbol{Y}^\text{T}-\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T})(\boldsymbol{Y}-\boldsymbol{X}\hat{\boldsymbol{\beta}}) \\ &=& \boldsymbol{Y}^\text{T}\boldsymbol{Y}-\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T}\boldsymbol{Y}-\boldsymbol{Y}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}+\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}} \\ &=& \boldsymbol{Y}^\text{T}\boldsymbol{Y}-2\boldsymbol{Y}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}+\hat{\boldsymbol{\beta}}^\text{T}\boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}} \end{array}$ (2-7)

然后對 D 求 $\hat{\boldsymbol{\beta}}$ 的偏導，并令其等于零咙轩，可得^[2]

$\begin{array}{rcl} \frac{\partial\boldsymbol{D}}{\partial\hat{\boldsymbol{\beta}}}=-2\boldsymbol{X}^\text{T}\boldsymbol{Y}+2\boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}&=&0\\ \Rightarrow \boldsymbol{X}^\text{T}\boldsymbol{X}\hat{\boldsymbol{\beta}}&=&\boldsymbol{X}^\text{T}\boldsymbol{Y}\end{array}$ (2-8)

式(2-8)即為法方程組的矩陣形式获讳。為了進一步證明極小值存在，可對D求 $\hat{\boldsymbol{\beta}}$ 的二階導數

$\frac{\partial^2\boldsymbol{D}}{\partial\hat{\boldsymbol{\beta}}^2}=2\boldsymbol{X}^\text{T}\boldsymbol{X}$ (2-9)

因此活喊，式(2-9)是正定陣丐膝，故該問題的最優(yōu)解存在。

接下來钾菊，分情況討論法方程組的解帅矗。
1）Rank( X )＝k＋1，即 X 列滿秩结缚。
令 C ＝ X^TX ，顯然软棺，Rank(C)＝k＋1红竭，且C是k階正定陣。因此喘落，C^-1存在茵宪，同時，可知式(2-9)正定瘦棋，故極值點存在稀火，得

$\hat{\boldsymbol{\beta}}=[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{Y}$ (2-10)

2）Rank( X )＝ r，0＜r＜k＋1赌朋。
此時凰狞，方程組并不符合假設3的條件篇裁，即存在完全共線的變量，但可以證明赡若，方程組是相容的^[3]达布，因此，方程組有唯一的極小范數解^[4]逾冬，為

$\hat{\boldsymbol{\beta}}=[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{+}\boldsymbol{X}]^\text{T}\boldsymbol{Y}$ (2-11)

這種情況下黍聂，估計量將無法保證是BLUE，此種情況已超出本文討論的范圍身腻，從略产还。此時，可將共線的變量合并后嘀趟，重新計算脐区。
由此也可看出，若分析場景中無法保證所選變量完全共線時去件，建議先通過其他相關性分析算法坡椒，對數據進行一次前期的清洗倔叼，選擇那些對因變量影響大宫莱、關聯(lián)性強的自變量丈攒，剔除那些關聯(lián)性弱或與其他自變量存在共線（或近似共線）的自變量授霸。

綜上所述，式(2-10)給出了計算最小二乘估計量的一般公式碘耳，根據此式显设，便可實現具體的最小二乘估計算法捕捂。

2.4斗搞、最小二乘估計量的統(tǒng)計特性

1）線性性
式(2-10)已表明僻焚， $\hat{\boldsymbol{\beta}}$ 是 Y 的線性組合。

2）無偏性
若證明 $\text{E}(\hat{\boldsymbol{\beta}})=\boldsymbol{\beta}$ 隙弛。
證明：

$\begin{array}{rcl}\text{E}(\hat{\boldsymbol{\beta}}) &=& \text{E}([(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{Y})\\&=&\text{E}([(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}(\boldsymbol{X}\hat{\boldsymbol{\beta}}+\boldsymbol{\varepsilon}))\\&=&\text{E}(\hat{\boldsymbol{\beta}}+[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{\varepsilon})\\ &=& \text{E}(\boldsymbol{\beta}) \\ &=& \boldsymbol{\beta}\end{array}$ (2-12)

證畢绩蜻。█

3）有效性
有效性是個定性的特性室埋，即估計量的方差越小姚淆，則說明此估計量越有效腌逢。而G-M定理已指出，最小二乘估計量是最優(yōu)的佳鳖，即 $\text{Var}(\hat{\boldsymbol{\beta}})$ 是所有估計量中最小的∠捣裕現在簡單證明這一結論妒蔚。
證明：
首先求 $\text{Var}(\hat{\boldsymbol{\beta}})$ ^[5]肴盏。

$\begin{array}{rcl}\text{Var}(\hat{\boldsymbol{\beta}}) &=& \text{E}([(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta})(\hat{\boldsymbol{\beta}}-\boldsymbol{\beta})^\text{T}]\\&=&\text{E}([(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\text{T}\boldsymbol{X}(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1})\\&=& [(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}\text{E}(\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\text{T})\boldsymbol{X}(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\\ &=& [(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}(\sigma^2\boldsymbol{I})\boldsymbol{X}(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1} \\ &=& \sigma^2(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\end{array}$ (2-13)

再證明式(2-13)是最小的贞绵。
令 $\tilde{\boldsymbol{\beta}}$ 是無偏的恍飘，則

$\text{E}(\boldsymbol{CY})=\text{E}(\boldsymbol{CX\beta}+\boldsymbol{C\varepsilon})=\boldsymbol{\beta}$

這表明常侣，CX＝I弹渔。把 C 代入到式(2-10)、式(2-13)中可得舞肆，

$\text{Var}(\tilde{\boldsymbol{\beta}})=\sigma^2\boldsymbol{CC}^\text{T}$

令 $\boldsymbol{D}=\boldsymbol{C}-[(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X}]^\text{T}$ 椿胯，則

$\text{Var}(\tilde{\boldsymbol{\beta}})=\sigma^2\boldsymbol{DD}^\text{T}+\sigma^2[(\boldsymbol{D}+(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X})^\text{T})(\boldsymbol{D}+(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}\boldsymbol{X})^\text{T})^\text{T}]$

由于 $\boldsymbol{DX}=\boldsymbol{\Theta}$ ，所以

$\text{Var}(\tilde{\boldsymbol{\beta}})=\sigma^2\boldsymbol{DD}^\text{T}+\sigma^2(\boldsymbol{X}^\text{T}\boldsymbol{X})^{-1}=\text{Var}(\hat{\boldsymbol{\beta}})+\sigma^2\boldsymbol{DD}^T$ (2-14)

因為DD^T是非負定矩陣^[6]前方，所以 $\text{Var}(\hat{\boldsymbol{\beta}})$ 是最小的惠险。█

綜上所述班巩，最小二乘估計量是最小方差線性無偏估計量（BLUE）抱慌。

2.5眨猎、OLS推廣

實際上宵呛，我們經常使用OLS對非線性系統(tǒng)進行回歸分析宝穗。通常逮矛，有以下幾種情形：多項式回歸、對數回歸鲸伴、分段回歸等汞窗。
1）多項式回歸
令x＝(x⁰, x¹, x², …, x^k)^T仲吏，代入方程組進行求解裹唆。
但需注意，多項式的階數不能過高劳坑，一般地距芬，當k>3時蔑穴，擬合誤差會變得很大惧浴，擬合結果基本無法使用衷旅。因此柿顶，對于復雜曲線的擬合嘁锯，可將曲線進行分段家乘，然后采用低階多項式進行擬合計算。
2）對數回歸
令z＝lnx（或z＝e^x）耀找，代入方程組進行求解野芒。
采用對數回歸時狞悲，應注意x的定義域摇锋，若?x≤0乱投，則不能使用此方法顷编。另外媳纬，對數有“縮小”效果钮惠，指數有“放大”效果素挽。當x的范圍跨度過大（lg(max(x)－min(x))＞3）预明，且無明顯線性特征時撰糠，可選用對數變換；當x的變化幅度過兄继弧（(max(x)－min(x))/mean(x)?1）砚尽，且無明顯線性特征時尉辑，可選用指數變換隧魄。

3购啄、OLS實現

對于簡單的低階（k≤2）估計而言狮含，直接使用式(2-10)即可几迄。但注意到，式(2-10)中需要計算矩陣 X^TX 的逆木羹，而隨著 k 的增長坑填，這一步將變得十分低效脐瑰，所以苍在，應考慮使用某種方法避開矩陣求逆忌穿，于是我們想到了矩陣的QR分解掠剑。

3.1朴译、QR分解

定理3.1 　設 $\boldsymbol{A}\in\mathbb{R}_n^{n\times n}$ 眠寿，^[7]則存在正交陣^[8] Q 及正線上三角陣 R 盯拱，使得 A ＝ QR 狡逢。
證明：因為 A 是滿秩方陣奢浑，將 A 寫成列分塊形式： A ＝(a₁, a₂, …, a_n)雀彼，則a₁, a₂, …, a_n 均為 $\frak{B}=\lbrace \boldsymbol{u}_1, \boldsymbol{u}_2,..., \boldsymbol{u}_n\rbrace$ 徊哑，且 (a₁, a₂, …, a_n)＝(u₁, u₂, …, u_n)R 成立莺丑，其中 R 為正線上三角陣窒盐。
但由于 $\small a_i\in\mathbb{R}^n,u_i\in\mathbb{R}^n$ 蟹漓，故有

$\boldsymbol{A} ＝ \boldsymbol{QR}$ (3-1)

注意 Q＝( u ₁, u ₂, …, u _n )的各列標準正交葡粒，故 Q^TQ＝I 嗽交，即 Q 為一個正交陣夫壁。█
Q有以下良好的性質：
1）Q 是非奇異的盒让，且det(Q)＝±1邑茄；
2）Q ^-1＝ Q ^T肺缕；
3）正交陣的乘積仍是正交陣同木；
4）Q 的特征值的模為1泉手。
定理3.2（Householder變換）　對于任意 $\boldsymbol{H}=\boldsymbol{I}_n-2\boldsymbol{uu}^\text{T}$ 斩萌，使得

$\boldsymbol{Hx}=a\boldsymbol{e}$ (3-2)

成立颊郎。其中 $\boldsymbol{e}$ 為n維自然基^[9]姆吭， $a$ 是實數内狸，且 $|a|=‖\boldsymbol{u}‖_2$ 昆淡。
計算可得昂灵，

$\boldsymbol{u}=\frac{\boldsymbol{x}-a\boldsymbol{e}}{\left\|\boldsymbol{x}-a\boldsymbol{e}\right\|_2}\ \ (|a|=\left\|\boldsymbol{x}\right\|_2)$ (3-3)

定理3.2說明眨补，任意列向量被Householder陣左乘后撑螺，可變換成與某自然基的共線向量甘晤。
推論1（Householder變換法） 對于滿秩方陣 A 安皱，將其按列分塊得 A ＝(α₁, α₂, …, α_n)酌伊，取

$\boldsymbol{u}_1=\frac{\boldsymbol{\alpha}_1-a_1\boldsymbol{e}_1}{\left\|\boldsymbol{\alpha}_1-a_1\boldsymbol{e}_1\right\|_2}\ \ (|a_1|=\left\|\boldsymbol{\alpha}_1\right\|_2)$ (3-4)

則 $\boldsymbol{H}_1=\boldsymbol{I}-2\boldsymbol{u}_1\boldsymbol{u}_1^\text{T}$ 虹脯，有

$\boldsymbol{H}_1\boldsymbol{A}=(\boldsymbol{H}_1\boldsymbol{\alpha}_1+\boldsymbol{H}_1\boldsymbol{\alpha}_2+\ldots+\boldsymbol{H}_1\boldsymbol{\alpha}_n)=\left[\begin{matrix} a_1 & * & \ldots &* \\0& \\ \vdots & &\boldsymbol{B}_1 & \\0& \end{matrix} \right]$ (3-5)

將矩陣 $\boldsymbol{B}_1\in\mathbb{R}^{(n-1)\times (n-1)}$ 按列分塊循集，得 B₁ ＝(β₂, …, β_n)咒彤，取

$\boldsymbol{u}_2=\frac{\boldsymbol{\beta}_2-a_2\boldsymbol{e}_2}{\left\|\boldsymbol{\beta}_2-a_2\boldsymbol{e}_2\right\|_2}\ \ (|a_2|=\left\|\boldsymbol{\beta}_2\right\|_2)$ (3-6)

$\tilde{\boldsymbol{H}}_2=\boldsymbol{I}-2\boldsymbol{u}_2\boldsymbol{u}_2^\text{T}\ \ \ \boldsymbol{H}_2=\left[\begin{matrix}1&\boldsymbol{\theta}^\text{T}\\\boldsymbol{\theta}&\tilde{\boldsymbol{H}}_2\end{matrix}\right]$ (3-7)

則

$\boldsymbol{H}_2\boldsymbol{H}_1\boldsymbol{A}=\left[\begin{matrix} a_1 & * & \ldots& \ldots &* \\0&a_2&*&\ldots&* \\ \vdots &0 \\ \vdots&\vdots&&\boldsymbol{C}_2 & \\0&0 \end{matrix} \right]$ (3-8)

依次進行下去，得到第n－1階的Householder陣H_n－1鞋屈，得

$\boldsymbol{H}_{n-1}\ldots\boldsymbol{H}_2\boldsymbol{H}_1\boldsymbol{A}=\left[\begin{matrix} a_1 & * & \ldots &* \\&a_2&\ldots&* \\ &&\ddots&\vdots\\&&&a_n \end{matrix} \right]=\boldsymbol{R}$ (3-9)

由于H^-1＝H厂庇，所以令Q＝H₁H₂…H_n-1权旷，則A＝QR炼杖。█

由此可見，使用推論1提供的方法罚缕，對于任意滿秩方陣邮弹，只需經過向量加減腌乡、矩陣相乘運算和基本四則運算就能實現QR分解与纽。

現將X進行QR分解急迂，即令X＝QR僚碎，則式(2-8)得

$\begin{aligned}\ \ \ \ \boldsymbol{X}^\text{T}\boldsymbol{X}\boldsymbol{\hat{\beta}}&=\boldsymbol{X}^\text{T}\boldsymbol{Y}\\ \Leftrightarrow (\boldsymbol{QR})^\text{T}\boldsymbol{QR}\boldsymbol{\hat{\beta}}&=(\boldsymbol{QR})^\text{T}\boldsymbol{Y}\\ \Leftrightarrow \boldsymbol{R}^\text{T}\boldsymbol{Q}^\text{T}\boldsymbol{QR}\boldsymbol{\hat{\beta}}&=\boldsymbol{R}^\text{T}\boldsymbol{Q}^\text{T}\boldsymbol{Y} \\ \Leftrightarrow \boldsymbol{R}^\text{T}\boldsymbol{R}\boldsymbol{\hat{\beta}}&=\boldsymbol{R}^\text{T}\boldsymbol{Q}^\text{T}\boldsymbol{Y} \\ \Leftrightarrow \boldsymbol{R}\boldsymbol{\hat{\beta}}&=\boldsymbol{Q}^\text{T}\boldsymbol{Y}\end{aligned}$ (3-10)

求解方程組(3-10)顯然不需要求 R^-1 了卷中，這是因為 R 是上三角陣仓坞，所以只需先計算出 β_n 无埃，然后將其代入到第n－1個方程中嫉称，計算出β_n－1；再將β_n震捣、β_n－1代入到第n－2個方程中蒿赢，計算出β_n-2壹若，以此類推店展，直至求出所有β為止赂蕴。

綜上所述概说，借助QR分解的方法席怪，可以避免矩陣求逆纤控，有效地降低了最小二乘算法的復雜度刻撒。

3.2声怔、Apache Commons Math

Apache公共包中有一個專門用于數學計算的工具包Math（org.apache.commons.math3）^[10]醋火，它提供了豐富且高效的數學分析算法芥驳，這其中就包含普通最小二乘法（org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression類）兆旬。
OLSMultipleLinearRegression類的求解分兩步：
1）按照推論1給出的方法，對X進行QR分解丽猬；
2）借助org.apache.commons.math3.linear.QRDecomposition類宿饱，求解方程組(3-10)。
此外脚祟，還提供評價OLS效果的一些系數計算方法谬以，如SSR、SSTO愚铡、r²胡陪，調整的r²沥寥。

3.3、應用實例

基于Math包柠座，我們可以輕易地開發(fā)出應用最小二乘法進行數據分析的功能邑雅。下面就以對一個時間序列，進行二次多項式擬合為例妈经，說明OLSMultipleLinearRegression類的具體使用方法淮野。
例：設某時間序列為 S＝(1.6, 4.4, 9.4, 16.6, 25.5)^[11]，猜測其符合s(k)＝b₀＋b₁k＋b₂k²吹泡，k＝1, 2, 3, 4, 5≈栊牵現使用最小二乘法對 S 進行二次多項式擬合。
解：構造
$\boldsymbol{X}=\left[\begin{matrix}1&1&1\\1&2&4\\1&3&9\\1&4&16\\1&5&25\end{matrix}\right],\ \boldsymbol{Y}=\left[\begin{matrix}1.6\\4.4\\9.4\\16.6\\25.5\end{matrix}\right]$
然后傳給如下代碼：

 OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
 regression.newSampleData(Y, X);
 double[] b = regression.estimateRegressionParameters();
 // s=b[0]+b[1]*k+b[2]*k*k
 double r2 = regression.calculateAdjustedRSquared();

所得數組b爆哑，即為擬合系數洞难。█

4、廣義最小二乘法（GLS）簡介

在實際情況中揭朝，OLS成立的5個假設常常難以滿足队贱。當假設5（同方差性）不能滿足時色冀，即ε的方差隨輸入的變化而變化，那么之前對OLS的推導過程就不再成立了柱嫌。對此锋恬，可以通過一些線性變換，使得模型滿足同方差的假設编丘，即又轉換成OLS模型了与学。
具體地，設線性系統(tǒng) Y ＝ Xβ ＋ ε 嘉抓，若Var(ε)＝σ²Ω癣防，其中，Ω是正定對稱陣掌眠。
則令Ω^-1＝P^TP蕾盯，則Ω＝P^-1(P^T)^-1，那么蓝丙，