劉應興的部落格

線性模型，以矩陣表示,可寫成  y = Xβ + ε。以 b 表示 β 之最小平方估計，則　

b = (X'X)^(-1)(X'y), 其中 X' 是 X 的轉置.

迴歸之殘差平方和為

SSE = y'y - b'X'Xb = y'y - y'X(X'X)^(-1)X'y,  或  y'(I-M)y

這算式其中一個重點是計算 X'X 的反矩陣，可有多種不同解法，基本的有解 AB = I 的消去法，即：對擴充矩陣 [A : I] 做基本列運算把前半部變成 I，即 [A : I] ～ [I : A^(-1)]。另外，有對 A 做 LU 分解、QR 分解、Cholesky 分解、singular value 分解等方法，藉由分解 A 成為特殊矩陣乘積，如 A = LU 而 AX = I 或一般 AX = B 變成 LUX = B，而解 LY = B 和 UX = Y 都是容易的。

不過, 在最小平方問題上, 卻另有一所謂 sweep 運算法, 將矩陣

[ X'X  X'y ] —————→ [ (X'X)^(-1)       b         ]
[ y'X  y'y ]  變成  [   -b'         y'y - b'X'Xb ]

結果包含估計的迴歸係數和誤差平方和。重要的是：如果 X 分成兩部分 X = [X1 : X2]，在計算過程將可看到資料配適子模型 y = X1．β1 + δ 的結果：(令 M1 = X1(X1'X1)^(-1)X1'，稱為在 C(X1) 的垂直投影算子)

　[ X1'X1  X1'X2  X1'y ]
　[ X2'X1  X2'X2  X2'y ]
　[  y'X1    y'X2    y'y  ]
  　    [     (X1'X1)^(-1)         (X1'X1)^(-1)X1'X2   (X1'X1)^(-1)X1'y ]
　  → [ -X2'X1(X1'X1)^(-1)     (X2'X2)-X2'M1X2       X2'y-X2'M1y      ]
  　    [ -y'X1(X1'X1)^(-1)          y'X2-y'M1X2            y'y-y'M1y         ]

為了簡化，我們用 A~F 來替代 X1'X1 等矩陣，並以 A^ 代表反矩陣，即 A^ = A^(-1), 則上列過程以至完結，是

[ A   B   D ]    [   A^        A^B          A^D    ]   令 C* =  (C-B'A^B)^ 
[ B'  C   E ] → [ -B'A^   C-B'A^B    E-B'A^D ]  ——————————→
[ D'  E'  F ]     [ -D'A^   E'-D'A^B   F-D'A^D ] 
    [    A^+A^BC*B'A^                   -A^BC*             A^D-A^BC*(E-B'A^D)            ]
    [          -C*B'A^                           C*                        C*(E-B'A^D)                   ]
    [ -D'A^+(E'-D'A^B)C*(B'A^)  -(E'-D'A^B)C*   (F-D'A^D)-(E'-D'A^B)C*(E-B'A^D) ]

在上列最後一個 3×3 區塊矩陣中，左上 2×2 區塊子矩陣其實是原區塊左上 2×2 區塊的反矩陣。也就是說，原模型

 y = X1．β1 + X2．β2 + ε 

對應的方陣經第一部分 sweep 運算得子模型 y = X1．β1 + δ 的迴歸係數估計及殘差平方和，又經後面部分的 sweep 運算，最後得全模型的迴歸係數估計及殘差平方和。

實際在做 sweep 運算是以 X'X 的主對角線元素當中樞 (pivot) 一個一個執行的，相當於對模型中的解釋變數及常數項一個一個加進去計算。首先是常數項，然後 X1，然後 X2，以此類推，至所有解釋變數都 sweep 為止。設欲做 sweep 運算的矩陣，以一般元素表示，是 [ a(ij) ]，以 a(kk) 為中樞做 sweep 運算後的矩陣是 [ b(ij) ]，則上述 sweep 運算的具體算式為：

b(kk) = 1/a(kk)
b(ik) = -a(ik)/a(kk),  i≠k
b(kj) = a(kj)/a(kk),   j≠k
b(ij) = a(ij)-a(ik)a(kj)/a(kk),  i≠k, j≠k

這算法並不是絕對的，例如可以 b(ik) = a(ik)/a(kk) 而 b(kj) = -a(kj)/a(kk)，則右邊得 （-b） 而下方是 b；或者，b(kk) = -1/a(kk) 而前列 b(ik) 及 b(kj) 均取正號，則最後結果是：

[ -(X'X)^(-1)      b       ]
[        b'         y'(I-M)y ]

也就是說：b(kk), b(ik), b(kj) 三類算式中恰一類取負號。這樣做的目的是在整個 sweep 運算過程中能夠一致。

個人感想：初見 sweep 運算法時，並不能理解其算法。它不像解 AX = B 對擴充矩陣 [ A : B ] 做基本列運算成等價方程式，所有計算過程看起來很是神秘，如上半部的 [ X'X X'y ] 變成 [ (X'X)^(-1) (X'X)^(-1)X'y ] 前者是取反矩陣，相當於以 (X'X)^(-1) 乘兩次；後者卻只左乘 (X'X)^(-1) 一次。雖然可以下式將 sweep 運算理解為列運算：

[      (X1'X1)^(-1)       0  0 ]     [ X1'X1  X1'X2  X1'y  I  0 ]
[ -X2'X1(X1'X1)^(-1)   I  0 ] ． [ X2'X1  X2' X2  X2'y  0  I ]
[ -y'X1(X1'X1)^(-1)    0  1 ]     [  y'X1    y'X2    y'y   0  0 ]
     [ I  (X1'X1)^(-1)X1'X2  (X1'X1)^(-1)X1'y     (X1'X1)^(-1)        0 ]
  = [ 0    X2'X2-X2'M1X2       X2'y-X2'M1y      -X2'X1(X1'X1)^(-1)  I ]
 　 [ 0     y'X2-y'M1X2            y'y-y'M1y         -y'X1(X1'X1)^(-1)   0 ]

最後是

[   (X'X)^(-1)      0 ] . [ X'X  X'y  I ] = [ I   (X'X)^(-1)X'y     (X'X)^(-1)    ]
[ -y'X(X'X)^(-1)  1 ]   [ y'X  y'y  0 ]    [ 0       y'y-y'My      -y'X(X'X)^(-1) ]

也就是將原來的平方和及交叉乘積和矩陣右邊填加含 I 及 0 的區塊，而後經列運算將左邊區塊變成 I 及 0，則後面兩區塊經行重排即是前面 sweep 運算結果。不過，這仍然像是在湊答案，並非真正的原理。究竟 sweep 運算的原理是怎麼回事，我仍不得其解。