第二章 文獻回顧
第五節 最小一乘法
國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
第五節 最小一乘法
最小一乘法(Least Absolute Deviation,LAD)，由 Boscovich 於 1755-1757 年提出，相較於Gauss 在 1809 年所提出的 LS 早了 40 年，但由於 LAD 的 目標函數為殘差絕對值之總和，其中的絕對值難以透過微分方式進行計算，
求解相當困難，導致LAD 在應用層面的發展停滯不前。但隨著科技發展、
演算法成熟，Charnes, Cooper W 和 Ferguson 等人於 1955 年，利用線性規 劃的方式進行求解，克服了LAD 目標函數難以微分的問題，將 LAD 應用 於管理問題中，而後LAD 才又開始被應用在不同領域中。
雖然LS 的演算簡單、速度快，使得 LS 被廣泛應用，但在觀測量含有 粗差時，會使LS 無法滿足高斯分布，導致成果受到影響(謝開貴等人，2002)；
而LAD 利用線性規劃的方式進行求解，雖然使得解算的複雜度提高，但此 種方式能夠克服觀測量含有粗差的問題(Nobakhti et. al,2009)，使成果不易 受粗差影響(王文峰，2006；Bektas,2010)，且具有較佳的穩健性。
許多不同的研究利用直線擬合來驗證"LAD 不易受粗差影響"的特性，
例如在實驗的數據中加入粗差，並對LAD 及 LS 的擬合成果進行評估，評 估結果發現LS 的迴歸直線並不通過任何一組數據點，而 LAD 則必通過其 中兩組數據點，越多的數據可以讓LAD 挑選出更佳的兩個通過點(李仲來，
1992)。當數據含有粗差時，LS 所計算出的成果與不含粗差的成果差異較 大，而透過LAD 所計算的兩個成果差異並不大，藉此展現 LAD 的穩健性 (謝開貴等人，2002)。
雖然LAD 具有穩健性，但缺乏統計上的指標，僅能以殘差進行比較，
並無法判斷可能為粗差的觀測量；在權迭代的計算方面，初值對於迭代結 果具有相當大的影響。因此，趙言等人於2016 年利用 LAD 計算的成果作 為權迭代法的初值，增加權迭代法的偵錯能力，其實驗結果顯示，該方法 除了能夠具有LAD 的穩健性外，亦能具有 LS 的統計指標，藉此彌補 LAD 無法對成果進行評估的缺陷。
針對 LAD 缺乏統計檢定量的問題，本研究提出一個新的概念”最小一 乘法的反求權矩陣”。在最小二乘的數學模型中，可分為函數模型及隨機模
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
精度的描述，而平差的目的即是求得未知參數的估計值，並估計各觀測量 的精度(李德仁，2002)。最小一乘法與最小二乘法的不同之處在於目標函數，
但所使用的函數模型是相同的，由前人文獻中可得知，當觀測量不含粗差 時，最小一乘法與最小二乘法所求得之期望值並無顯著性差異，但在觀測 量含有粗差時，最小一乘法的結果則較不易受到影響，而最小二乘法則須 配合適當的隨機模型才能降低粗差的影響。因此，本研究將最小一乘法視 為最小二乘法透過某種隨機模型所求出的期望值，透過最佳化演算法計算 LAD 的反求權矩陣，並利用此隨機模型所賦予的統計檢定量與其他方法進 行分析。
‧
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
其分母𝜎為標準差，為一固定值，因此分子部份必須為最小值，進而得出 (x − μ)2 = min ( 13 ) x 即為觀測量、μ即為最優估值、(x − μ)即為殘差 V，當殘差平方和最小時，
為出現的機率最高的狀況。因此，可將LS 的觀測方程式( 8 )進一步改寫為：
V = L − AX ( 14 ) 並以殘差平方之總和最小為目標函數，可對目標函數進一步推導
VTPV = min
VTPV = (L − AX)TP(L − AX)
= (LT− XTAT)P(L − AX)
= LTPL − LTPAX − XTATPL + XTATPAX ( 15 ) 而後對X 偏微分求VTPV極值可得
∂VTPV
∂X = −2ATPL + 2ATPAX = 0 ( 16 )
移項後可得
ATPAX = ATPL ( 17 ) 令式( 17 )中的
N = ATP𝐴 U = ATPA
N 為一n × n的方陣，可求得逆矩陣N−1，得式( 18 )：
X = N−1U ( 18 )
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
第二節 李德仁權迭代法
由文獻回顧可知，由於僅用改正數無法準確的判斷粗差的位置，若是 以此為權函數則有可能誤判觀測量。因此，將觀測量納入隨機模型，透過 權函數判斷觀測量精度後，降低可能為粗差的觀測量之權重，藉此降低該 觀測量的影響。
李德仁(1984)利用 Baarda 於 1968 年所提出的 Data Snooping 的概念，
將標準化殘差做為權迭代法中的權函數，而標準化殘差與多餘觀測分量(ri) 及單位權中誤差(σ0)有關。因此，在迭代的過程中σ0亦為影響的因素之一，
若只以最初的先驗中誤差做為判斷標準，容易有誤判的情形發生。因此，
應由前一次計算出的後驗中誤差(σ̂0)在權函數中進行取代，且在迭代初期應 給予較嚴格的限制，當迭代次數 n 大於 3 後再放寬限制，藉由統計檢定的 方式，對含有粗差的觀測量進行降降權，降低粗差對未知參數的影響，使 求解參數從有偏估值至無偏估值，其權函數如式( 7 )。
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
第三節 最佳化權矩陣
由LS 所得之X矩陣為：
X = (ATPA)−1ATPL ( 19) 將其代入式( 14 )可得
V = L − A(ATPA)−1ATPL ( 20 ) 當觀測量不含粗差時，權矩陣的變化並不會影響X 矩陣的期望值E(X̂)，
X̂為無偏估值；但當觀測量包含粗差時，不同的權矩陣則會影響估計參數的 求解。因此，各觀測量精度不同，應給予不同的權重，以降低含有粗差的 觀測量對最終求解成果的影響。
不同於權迭代法，最佳化權矩陣利用最佳化演算法進行計算，以標準 化殘差之總和最小為目標函數，求得權重於何種分布時為數據最集中的狀 態，而非透過統計檢定的方式在迭代中對觀測量進行降權。
以下將針對本研究所使用的最佳化演算法PSO 及優化後的 QPSO 進行 介紹。
一、 粒子群優化演算法
粒子群優化演算法(Particle Swarm Optimization, PSO)為 Kennedy &
Eberhart 於 1995 年提出基於群體智慧所進行的演算法。PSO 利用多個粒子 所組成的粒子群，在給定的搜尋空間(search space)中搜尋粒子的最佳位置，
並將其視為整體最佳解。在PSO 演算的過程中，每個粒子的位置都被視為 一組目標函數(objective function)的解，在迭代過程中不斷地更新粒子位置 及速度，藉此在空間中找到粒子最佳位置，並將其視為目標函數的最佳解。
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
PSO 目標函數為f(x)，且假設在 D 維的搜尋空間下，每一維度都有 M 個粒子組成立子群，每個粒子為Xi,1, Xi,2, … Xi,D, i = 1,2, … , M ;接著由式( 21 ) 更新粒子的速率及位置：
Vij(t + 1) = Vi,j(t) + c1∙ r1,i,j(t) ∙ (Pi,j(t) − Xi,j(t)) + c2∙ r2,i,j∙ (Gj(t) − Xi,j(t)) Xi,j(t + 1) = Vi,j(t + 1) + Xi,j(t) ( 21 ) 其中，Vij為粒子的速率；t 為迭代次數；Pi,j為粒子所搜尋到的最佳位置；
Gj為粒子群搜尋到的最好位置；r1,i,j與r2,i,j為兩個獨立且介於0~1 的隨機數；
c1、c2為個體與群體間互動的相關係數。每次迭代後，各粒子更新個體最佳 位置的判斷式為：
Pi,d(t + 1) = { Pi,d(t), if f[Xi,d(t)] ≥ f[Pi,d(t + 1)]
Xi,d(t + 1), if f[Xi,d(t + 1)] < f[Pi,d(t + 1)]
( 22 ) 而粒子群所找到的全體最佳位置的公式為：
g = arg min {f[Pi,d(t)]}, g ∈ {1,2, … , m} ( 23 ) Gd(t) = Pg,d(t) ( 24 ) PSO 流程如圖 4 所示。
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
初始化粒子群
迭帶次數 for t = 1 to n
計算各粒子的 當前適應值
粒子當前適應值 是否優於前一次結果 演算法開始
是
粒子更新坐標位置
否 粒子位置維持不變
更新群體最佳位置
利用速度及位置更新公式 更新粒子速度及位置
是否達到最大迭
否 帶次數 輸出成果
演算法結束
圖 4 PSO 流程
‧
量子行為粒子群演算法(Quantum-behaved Particle Swarm Optimization, QPSO)為孫俊等人於 2004 提出，為 PSO 的優化演算法，在 PSO 中引入量
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
其中，mbestd為 d 維空間中所有粒子的最佳位置的平均值；Pi,d為第 i 個粒子在 d 維空間中的最佳位置；𝐺𝑑為群體最佳位置；𝑝i,d為𝐺𝑑及Pi,d之間 的吸引子(local attractor)；𝜑為常態分布之下介於 0~1 之間的隨機數；u 則是 另一個常態分布之下介於 0~1 之間的隨機數；而α即為前述的收縮擴張係 數。
將 QPSO 最終的粒子最佳位置視為該目標函數的最佳解，即本研究的 最佳化權矩陣。
QPSO 的計算過程如下：
1. 隨機給予各粒子初始位置。
2. 利用上式計算各維度的平均最好位置𝐦𝐛𝐞𝐬𝐭𝐝。 3. 決定各粒子的當前最佳位置𝐏𝐢,𝐝。
4. 決定各維度的群體最佳位置𝑮𝒅。
5. 比較迭代後與迭代前一次之群體最佳位置，若是優於前一次迭代 值，則進行取代，否則維持原位置。
6. 計算介於𝑮𝒅與𝐏𝐢,𝐝間的隨機吸引子。
7. 利用步驟6 計算的吸引子，計算各粒子的更新位置。
8. 重複步驟2~7，直到迭代次數結束。
QPSO 流程如圖 5。
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
初始化粒子群
迭帶次數 for t = 1 to n
計算各粒子的 當前適應值
粒子當前適應值 是否優於前一次結果 演算法開始
是 粒子更新坐標位置
否 粒子位置維持不變
更新群體最佳位置
計算每個粒子的新位置
是否達到最大迭 否 帶次數
輸出成果
演算法結束 計算平均最佳位置
mbest
計算隨機吸引子
‧
‧
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
目標函數的係數矩陣（Cost Vector） c = [c1 c2 … cn]T 決策變數矩陣（Solution Vector） x = [x1 x2 … xn]T
右端項係數矩陣（Right-hand-side Vector） b = [b1 b2 … bn]T 約束條件的係數矩陣（Constraint Matrix） A = [
𝑎11 ⋯ 𝑎1𝑛
⋮ ⋱ ⋮ 𝑎𝑚1 ⋯ 𝑎𝑚𝑛] 二、線性規劃解
由目標函數及約束條件等線性條件組成可行域(F)，並於可行域中找出 最優解即為線性規劃問題。而線性規劃具有的特性為(方述誠，1993；林怡 君，2013)：
(1)線性規劃的可行解的集合，若非空即合則含有頂點。
(2)而各個頂點可能為最佳解。
因此線性規劃的解法可分成下列三個步驟(方述誠，1993)：
步驟一： 找出可行域(F)的各個頂點。
步驟二： 計算各個頂點所對應的目標函數值，即CT𝑋。
步驟三： 找出目標函數值最小的頂點，該頂點即為線性規劃的最優解。
單行法由George Dantzig 於 1947 年提出，為線性規劃中最常被使用、
且較有效率的方法之一，其構想為，先找尋可行域中的一個頂點，並比較 各個鄰近頂點的目標函數值，若有目標函數值更低的頂點，則往該頂點移 動，依照此規則依序檢視頂點，最終將找出目標函數值較鄰近頂點低或相 當的頂點(方述誠，1993)。圖 6 為單行法示意圖，在可行解中有多個頂點 (x1、x2…)，單行法由其中一個頂點(xk)開始搜索，沿著邊界移動到另一個 頂點(xk+1)，值到找到最佳解(x∗)為止。
圖 6 單行法示意圖(方述誠，1993)
‧
‧ 國
立 政 治 大 學
‧
N a tio na
l C h engchi U ni ve rs it y
而約制條件中，只需設計兩個部分，一個是觀測方程式中的等式L = AX，
第二個是待待求解矩陣中的各項元素皆大於等於0。因此：
Aeq = [ 𝐴 − 𝐴 − 𝐼 𝐼 ]、beq = 𝐿 lb = [
0
⋮ 0
]
將設計完成的矩陣依格式輸入，即可求得x 矩陣，再將X+減去X−、V+ 減去V−，即可求得未知參數X及改正數V。
四、結合最小一乘法與權迭代法
由於在迭代計算的過程中，給予的初值對於最終收斂結果具有相當大 的影響，不好的初值可能造成結果不如預期。因此，為了使迭代過程中的 權函數能夠更合理的反應隨機模型，趙言等人(2016)利用 LAD 解做為權迭 代法的初值，解決 LS 在含有粗差時容易因初值不正確而導致權迭代法成
由於在迭代計算的過程中，給予的初值對於最終收斂結果具有相當大 的影響，不好的初值可能造成結果不如預期。因此，為了使迭代過程中的 權函數能夠更合理的反應隨機模型，趙言等人(2016)利用 LAD 解做為權迭 代法的初值，解決 LS 在含有粗差時容易因初值不正確而導致權迭代法成