我在上文「對最小二乘法擬合直線的一些思考與改進」中介紹了通過距離的平方來擬合一條直線，這實際上也是一個PCA問題，但是我在結尾也說了，這種擬合方法的缺點是對離群噪點的幹擾異常明顯，稍不留神就把我們的目標給拉偏了，比如下面這張圖片，在工業機器視覺檢測中，為了測量某些產品的高精密尺寸，我們需要擬合出圖中斜邊邊緣的準確直線方程，但是由於產品的缺陷，在邊緣處往往有噪點的存在：
而如果我們不加以處理而直接擬合的話，最終的擬合結果是這樣的：
綠線是我們根據上篇文章推導出來的公式進行擬合出來的結果，這顯然不是我們所要的結果。由於噪點的幹擾，我們所需的邊緣直線方程被拉偏了。
這是由於在實際擬合過程中，離群較大的點它們都有很大的權重，而實際上每個點都應該有不同的權重，越偏離目標點它們的權重應該越小。因此我們迫切需要對每個點都定義一個權重，這就是今天我要介紹的加權最小二乘法。
在介紹這個算法之前，先回答一個問題，上篇文章中，有網友私信問我為什麼那個參數方程要選取較小的特徵值與特徵向量：
這個問題我在上篇文中只是提到了一下，最小距離的總和剛好是數據集矩陣的特徵值，但並沒有說明具體的原因，這裡給出一個更直觀的證明方法：
如上圖所示，我們定義一個直線的單位法向量
也即：
因此：
好了，我們再回到剛開始那個問題，我們究竟該對每個點如何定義權重值
但是我們究竟如何這個權重閾值
那麼我們該採用哪種方式進行擬合呢，上面說了，權重的取值跟點到直線的距離密切相關，因此我們依然採用點到直線距離的最小來擬合，這也是一個條件極值問題，同樣需要適配拉格朗日乘子：
下面給出實例代碼(c++版)：
bool FitLineWeight(vector<POINT> vec, double &a, double &b, double &c){ if (vec.size() < 2) { return false; } FitLine(vec, a, b, c);
for (size_t it = 0; it < times; it++) { if (vec.size() < 2) { return false; } double sum_x = 0.0f, sum_y = 0.0f; double avg_x = 0.0f, avg_y = 0.0f; vector<double> vecDistance; for (size_t i = 0; i < vec.size(); i++) { double d = abs(a * vec[i].x + b * vec[i].y + c); vecDistance.push_back(d); } sort(vecDistance.begin(), vecDistance.end());
double median; if (vecDistance.size() % 2 == 0) { median = (vecDistance[vecDistance.size() / 2 - 1] + vecDistance[vecDistance.size() / 2]) / 2; } else { median = vecDistance[vecDistance.size() / 2]; } double t = 2 * (median / 0.6745); double sum_ww = 0.0f; vector<double> vecW; for (int i = 0; i < vec.size(); i++) { double w; double d = abs(a * vec[i].x + b * vec[i].y + c); if (d <= t) { w = 1; } else { w = t / d; }
vecW.push_back(w); sum_x += w * w * vec[i].x; sum_y += w * w * vec[i].y; sum_ww += w * w; } avg_x = sum_x / sum_ww; avg_y = sum_y / sum_ww;
double L_xx = 0.0f; double L_yy = 0.0f; double L_xy = 0.0f; for (int i = 0; i < vec.size(); i++) { L_xx += pow(vecW[i], 2) * (vec[i].x - avg_x) * (vec[i].x - avg_x); L_xy += pow(vecW[i], 2) * (vec[i].x - avg_x) * (vec[i].y - avg_y); L_yy += pow(vecW[i], 2) * (vec[i].y - avg_y) * (vec[i].y - avg_y); }
double lamd1, lamd2; double lamd; double m, n; m = L_xx + L_yy; n = L_xx * L_yy - L_xy * L_xy; lamd1 = (m + sqrt(m * m - 4 * n)) / 2; lamd2 = (m - sqrt(m * m - 4 * n)) / 2; lamd = lamd1 < lamd2 ? lamd1 : lamd2;
double d = sqrt((L_xx - lamd) * (L_xx - lamd) + L_xy * L_xy); if (abs(d) < 1e-6) { a = 1; b = 0; c = -a * avg_x - b * avg_y; } else { if (lamd >= L_xx) { a = L_xy / d; b = (lamd - L_xx) / d; c = -a * avg_x - b * avg_y; } else { a = -L_xy / d; b = (L_xx - lamd) / d; c = -a * avg_x - b * avg_y; } } }}
按照上面的代碼，分別調整代碼裡面的循環變量
不迭代結果(也就是沒有權重分析，直接擬合，讓變量
讓
迭代兩次又更好了：
迭代三次的結果：
可見，只經過了三次迭代，我們擬合的直線就收斂於我們想要的理想結果了。而事實也證明，基本上只要我們迭代三到五次，就能達到很好的效果。
說到這裡，本文就算結束了，等我有空會再通過這種方法驗證一下其他集合圖形的擬合結果，例如圓或者橢圓，到時候再給大家回復。如果你想要本文的原始碼，請掃描下面的圖片，在我公眾號內回復「加權最小二乘」六個字即可，我會將原始碼都發送給你，方便你下去研究學習。