lm opencv 算法_LM算法
最小二乘法的概念
最小二乘法要關心的是對應的cost function是線性還是非線性函數，不同的方法計算效率如何，要不要求逆，矩陣的維數
一般都是過約束，方程式的數目多于未知的參數數目。
最小二乘法的目標：求誤差的最小平方和，根據cost function的對應有兩種：線性和非線性(取決于對應的殘差(residual)是線性的還是非線性的)。
線性最小二乘的解是closed-form solution 即 \(x = (A^TA)^{-1}A^Tb\)
\(S(x) = ||b -Ax ||^2 = (b -Ax)^T(b -Ax) = b^Tb - b^TAx - x^TA^Tb + x^TA^TAx\)
Note that :\((b^TAx)^T = x^TA^Tb\) 的維數是 1*1(y的列數目).所以 \(b^TAx = x^TA^Tb\)
\(S(x) = b^Tb - 2b^TAx + x^TA^TAx\)
對\(x\)求導
\(-2A^Tb + 2 (A^TA)x = 0\)
\(x = (A^TA)^{-1}A^Tb\)
而非線性最小二乘沒有closed-form，通常用迭代法求解。在每次迭代的過程使用一個線性化的方程代替計算。
有牛頓法，牛頓高斯法，LM, 其實可以分為trust region 和 linear line search
非線性最小二乘法的方法有
迭代法，即在每一步update未知量逐漸逼近解，cost function在下降，可以用于各種各樣的問題。
梯度下降(最速下降法)是迭代法的一種，可以用于求解最小二乘問題
\(f(x + \alpha \overrightarrow d) = f(x_0) + \alpha f'(x_0)\overrightarrow d + \text {二階以上無窮小}\)
上面的\(x_0\)是泰勒展開式的點， \(\alpha\)是 步長(一個實數) , $ \overrightarrow d$ 單位方向(一個向量)，即 \(|\overrightarrow d| = 1\)
顯然當\(d\)的方向和\(f'(x_0)\)方向相反的時候，\(cos180 =-1\),整體取到最小值。這就是為什么取的是梯度的反方向原因。
當然上面的步長一般也是可以求得，怎么求步長也是一門學問。
下面的都是從牛頓法引申出來的,記住牛頓法求得是穩定點\(f'(x) = 0\),導數為0的不一定是最小值,梯度下降法求得是局部最小值,從計算上看
牛頓法 泰勒展開式到二階，
\[f(x_{t+1}) = f(x_t) + g(x_{t+1} - x_t) + \frac{1}{2} (x_{t+1} -x_t)^TH(x_{t+1} -x_t)
\]
求導就有，\(\frac{\partial f}{\partial x_t} = g + H(x_{t+1} -x_t)\)，讓他為0，就有了牛頓公式
\[x_{t+1} = x_t - H^{-1}g
\]
要是H是正定的，上面的就是凸函數，也就一定有了最小值。可惜H不一定是正定的，這就引導出了下面的方法
高斯-牛頓法
是另一種經常用于求解非線性最小二乘的迭代法(一定程度上可視為標準非線性最小二乘求解方法)。
我們優化的cost function是 $$min \sum _i r_i(x)^2$$
根據上面的牛頓法：
\[x_{t+1} = x_t - H^{-1}g
\]
梯度的表示
\[g_j=2 \sum _i r_i \frac{\partial r_i}{\partial x_j}
\]
Hessian矩陣的表示
\[H_{jk} = 2\sum _i (\frac{\partial r_i}{\partial x_j}\frac{\partial r_i}{\partial x_k} + r_i\frac{\partial ^2 r_i}{\partial x_j \partial x_k})
\]
要是把上面公式的最后一項去掉，至少是半正定了，而且不用計算Hessain矩陣了，這就是牛頓高斯法
\[H_{jk} \approx 2 \sum _i J_{ij}J_{ik} \quad With \quad J_{ij} = \frac{\partial r_i}{\partial x_j}
\]
在什么情況下，去掉最后一項比較好，
residual(\(r_i\))小的時候
或者接近linear ,這樣一階微分是常數，二階微分就是0
上面兩種情況下，第二項都很小,公式如下
\[x_{t+1} = x_t - (J^TJ)^{-1}J^Tr
\]
Levenberg-Marquardt
的迭代法用于求解非線性最小二乘問題，就結合了梯度下降和高斯-牛頓法。
\[x_{t+1} = x_t - (H + \lambda I_n)^{-1}g
\]
\[x_{t+1} = x_t - (J^TJ + \lambda I_n)^{-1}J^Tr
\]
總結
所以如果把最小二乘看做是優化問題的話，那么梯度下降是求解方法的一種，\(x=(A^TA)^{-1}A^Tb\)是求解線性最小二乘的一種，高斯-牛頓法和Levenberg-Marquardt則能用于求解非線性最小二乘。
LM算法相對于高斯牛頓算法和梯度下降的優缺點
首先梯度下降法和高斯牛頓法都是最優化方法。其區別之處在于，
梯度下降法在尋找目標函數極小值時，是沿著反梯度方向進行尋找的。梯度的定義就是指向標量場增長最快的方向，在尋找極小值時，先隨便定初始點(x0，y0)然后進行迭代不斷尋找直到梯度的模達到預設的要求。但是梯度下降法的缺點之處在于：在遠離極小值的地方下降很快，而在靠近極小值的地方下降很慢，靠近的時候可能成zig-zag下降。
而高斯牛頓法是一種非線性最小二乘最優化方法。其利用了目標函數的泰勒展開式把非線性函數的最小二乘化問題化為每次迭代的線性函數的最小二乘化問題。高斯牛頓法的缺點在于：若初始點距離極小值點過遠，迭代步長過大會導致迭代下一代的函數值不一定小于上一代的函數值。
LM算法在高斯牛頓法中加入了因子μ，當μ大時相當于梯度下降法，μ小時相當于高斯牛頓法。在使用Levenberg-Marquart時，先設置一個比較小的μ值，當發現目標函數反而增大時，將μ增大使用梯度下降法快速尋找，然后再將μ減小使用牛頓法進行尋找。
The Gauss–Newton algorithm is used to solve non-linear least squares problems.
It is a modification of Newton's method for finding a minimum of a function.
Unlike Newton's method, the Gauss–Newton algorithm can only be used to minimize a sum
of squared function values, but it has the advantage that second derivatives, which
can be challenging to compute, are not required.
However, as for many fitting algorithms, the LMA finds only a local minimum,
which is not necessarily the global minimum.
% 計算函數f的雅克比矩陣
syms a b y x real;
f=a*cos(b*x) + b*sin(a*x)
Jsym=jacobian(f,[a b])
data_1=[ 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2 ];
obs_1=[102.225 ,99.815,-21.585,-35.099, 2.523,-38.865,-39.020, 89.147, 125.249,-63.405, -183.606, -11.287,197.627, 98.355, -131.977, -129.887, 52.596, 101.193,5.412, -20.805, 6.549, -40.176, -71.425, 57.366, 153.032,5.301, -183.830, -84.612, 159.602, 155.021, -73.318, -146.955];
% 2. LM算法
% 初始猜測初始點
a0=100; b0=100;
y_init = a0*cos(b0*data_1) + b0*sin(a0*data_1);
% 數據個數
Ndata=length(obs_1);
% 參數維數
Nparams=2;
% 迭代最大次數
n_iters=60;
% LM算法的阻尼系數初值
lamda=0.1;
%LM算法的精度
ep=100
% step1: 變量賦值
updateJ=1;
a_est=a0;
b_est=b0;
% step2: 迭代
for it=1:n_iters
if updateJ==1
% 根據當前估計值，計算雅克比矩陣
J=zeros(Ndata,Nparams);
for i=1:length(data_1)
J(i,:)=[cos(b_est*data_1(i))+data_1(i)*b_est*cos(a_est*data_1(i)) -sin(b_est*data_1(i))*a_est*data_1(i)+sin(a_est*data_1(i)) ];
end
% 根據當前參數，得到函數值
y_est = a_est*cos(b_est*data_1) + b_est*sin(a_est*data_1);
% 計算誤差
d=obs_1-y_est;
% 計算(擬)海塞矩陣
H=J'*J;
% 若是第一次迭代，計算誤差
if it==1
e=dot(d,d);
end
end
% 根據阻尼系數lamda混合得到H矩陣
H_lm=H+(lamda*eye(Nparams,Nparams));
% 計算步長dp，并根據步長計算新的可能的\參數估計值
dp=inv(H_lm)*(J'*d(:));
%求誤差大小
g = J'*d(:);
a_lm=a_est+dp(1);
b_lm=b_est+dp(2);
% 計算新的可能估計值對應的y和計算殘差e
y_est_lm = a_lm*cos(b_lm*data_1) + b_lm*sin(a_lm*data_1);
d_lm=obs_1-y_est_lm;
e_lm=dot(d_lm,d_lm);
% 根據誤差，決定如何更新參數和阻尼系數
if e_lm
if e_lm
break
else
lamda=lamda/5;
a_est=a_lm;
b_est=b_lm;
e=e_lm;
disp(e);
updateJ=1;
end
else
updateJ=0;
lamda=lamda*5;
end
end
%顯示優化的結果
a_est
b_est
plot(data_1,obs_1,'r')
hold on
plot(data_1,a_est*cos(b_est*data_1) + b_est*sin(a_est*data_1),'g')
#pragma once
#include
#include
#include
using namespace std;
using namespace cv;
const int MAXTIME = 50;
#pragma comment(lib,"opencv_core249d.lib")
Mat cvSinMat(Mat a)
{
int rows = a.rows;
int cols = a.cols;
Mat out(rows, cols, CV_64F);
for (int i = 0; i < rows; i++)
{
out.at(i, 0) = sin(a.at(i, 0));
}
return out;
}
Mat cvCosMat(Mat a)
{
int rows = a.rows;
int cols = a.cols;
Mat out(rows, cols, CV_64F);
for (int i = 0; i < rows; i++)
{
out.at(i, 0) = cos(a.at(i, 0));
}
return out;
}
Mat jacobin(const Mat& pk, const Mat& x) // pk= [a,b] a*cos(b*x) + b*sin(a*x)
{
Mat_ J(x.rows, pk.rows), Sa, Sb ,Ca, Cb,da, db;
Sa = cvSinMat(pk.at(0)*x);
Sb = cvSinMat(pk.at(1)*x);
Ca = cvCosMat(pk.at(0)*x);
Cb = cvCosMat(pk.at(1)*x);
da = Cb + x.mul(pk.at(1)*Ca);
db = Sa - x.mul(pk.at(0)*Sb);
//cout << "da= " << da << endl;
da.copyTo(J(Rect(0, 0, 1, J.rows)));
db.copyTo(J(Rect(1,0, 1, J.rows)));
return J;
}
Mat yEstimate(const Mat& pk, const Mat& x)
{
Mat_ Y(x.rows, x.cols),Cb,Sa;
Sa = cvSinMat(pk.at(0)*x);
Cb = cvCosMat(pk.at(1)*x);
Y = pk.at(0)*Cb + pk.at(1)*Sa;
return Y;
}
void LM(double* p0, int pN, double* x, int xN, double* y, double lamda, double step, double ep = 0.000001)
{
int iters = 0;
int updateJ = 1;
double ek = 0.0, ekk = 0.0;//估計誤差
Mat_ xM(xN, 1, x), yM(xN, 1, y), pM(pN, 1, p0), JM, yEM, yEMM, dM, gM, dMM, dpM;//至少需要JM,gM,dpM,pM
for (; iters < MAXTIME; iters++)
{
if (updateJ == 1)
{
JM = jacobin(pM, xM); //雅克比矩陣
//outData(fs, JM, "J.xml");
yEM = yEstimate(pM, xM); //f(β)
dM = yM - yEM; // y-f(β)
gM = JM.t()*dM; //
if (iters == 0)
ek = dM.dot(dM); //第一次 直接||r||^2
}
Mat_ NM = JM.t()*JM + lamda*(Mat::eye(pN, pN, CV_64F)); //J(T)J + lambda*I
if (solve(NM, gM, dpM)) //
{
Mat_ pMM = pM + dpM; //更新最小值
yEMM = yEstimate(pMM, xM);
dMM = yM - yEMM;
ekk = dMM.dot(dMM);
if (ekk < ek)//成功則更新向量與估計誤差
{
printf("the %d iterator ,ekk=%lf \n", iters,ekk);
if (dpM.dot(dpM) < ep)
{
printf("Final result is :\n");
printf("精度：0.000001\n");
printf("a =%lf , b=%lf\n", pMM.at(0),pMM.at(1));
return;
}
else
{
pM = pMM;
ek = ekk;
lamda = lamda / step;
updateJ = 1;
continue;
}
}
else //if an iteration gives insufficient reduction in the residual, λ(lamda) can be increased
{
printf("the %d iterator \n", iters);
lamda = lamda*step;
updateJ = 0;
}
}
else
{
printf("the solve invertx matrix error\n");
}
}
}
#include "LM.h"
int main()
{
double data[] = { 0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0,
3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2 };
double obs[] = { 102.225 ,99.815,-21.585,-35.099, 2.523,-38.865,
-39.020, 89.147, 125.249,-63.405, -183.606, -11.287,
197.627, 98.355, -131.977, -129.887, 52.596, 101.193,
5.412, -20.805, 6.549, -40.176, -71.425, 57.366, 153.032,
5.301, -183.830, -84.612, 159.602, 155.021, -73.318, -146.955 };
//初始點
double p0[] = { 100, 100 };
LM(p0, 2, data, 32, obs, 0.1, 10);
system("pause");
}
LM算法與非線性最小二乘問題
摘錄的一篇有關求解非線性最小二乘問題的算法--LM算法的文章,當中也加入了一些我個人在求解高精度最小二乘問題時候的一些感觸: LM算法,全稱為Levenberg-Marquard算法,它可用于解決非線 ...
Levenberg-Marquardt迭代(LM算法)-改進Guass-Newton法
1.前言 ?????????????????????????????? a.對于工程問題,一般描述為:從一些測量值(觀測量)x 中估計參數 p?即x = f(p), ??? ...
梯度下降法、牛頓法、高斯牛頓法、LM最優化算法
1.梯度下降法 2.牛頓法 3.高斯牛頓法 4.LM算法
Levenberg-Marquardt優化算法以及基于LM的BP-ANN
一.LM最優化算法 ? ? 最優化是尋找使得目標函數有最大或最小值的的參數向量.根據求導數的方法,可分為2大類.(1)若f具有解析函數形式,知道x后求導數速度快.(2)使用數值差分來求導數.根據使用模 ...
LM擬合算法
一.? Levenberg-Marquardt算法 (1)y=a*e.^(-b*x)形式擬合 clear all % 計算函數f的雅克比矩陣,是解析式 syms a b y x real; f=a*e ...
Levenberg-Marquardt算法基礎知識
Levenberg-Marquardt算法基礎知識 (2013-01-07 16:56:17) 轉載▼ ? 什么是最優化?Levenberg-Marquardt算法是最優化算法中的一種.最優化是尋找使 ...
相機標定：關于用Levenberg-Marquardt算法在相機標定中應用
LM算法在相機標定的應用共有三處. (1)單目標定或雙目標定中,在內參固定的情況下,計算最佳外參.OpenCV中對應的函數為findExtrinsicCameraParams2. (2)單目標定中,在 ...
點云匹配和ICP算法概述
Iterative Closest Point (ICP) [1][2][3] is an algorithm employed to minimize the difference between ...
隨機推薦
Hibernate中事務聲明
Hibernate中JDBC事務聲明,在Hibernate配置文件中加入如下代碼,不做聲明Hibernate默認就是JDBC事務. 一個JDBC 不能跨越多個數據庫. Hibernate中JTA事務聲 ...
Redis基礎知識之————空間換時間的查詢案例
空間與時間 空間換時間是在數據庫中經常出現的術語,簡單說就是把查詢需要的條件進行索引的存儲,然后查詢時為O(1)的時間復雜度來快速獲取數據,從而達到了使用空間存儲來換快速的時間響應!對于redis這個 ...
Basic Printing Architecture
https://blogs.technet.microsoft.com/askperf/2007/06/19/basic-printing-architecture/ Printer sharing, ...
Android 文字繪制(DrawText)技術總結
這里的繪制文字不是直接調用TextView.setText(String content)去展示文字內容.而是在View上面通過?canvas.drawText(text, x, y,textPain ...
iOS學習——如何在mac上獲取開發使用的模擬器的資源以及模擬器中每個應用的應用沙盒
如題,本文主要研究如何在mac上獲取開發使用的模擬器的資源以及模擬器中每個應用的應用沙盒.做過安卓開發的小伙伴肯定很方便就能像打開資源管理器一樣查看我們寫到手機本地或應用中的各種資源,但是在iOS開發 ...
洛谷P3835 【模板】可持久化平衡樹
題目背景 本題為題目?普通平衡樹?的可持久化加強版. 數據已經經過強化 感謝@Kelin 提供的一組hack數據 題目描述 您需要寫一種數據結構(可參考題目標題),來維護一些數,其中需要提供以下操作( ...
Stream初步應用
一.什么是stream Stream(流)是一個來自數據源的元素隊列并支持聚合操作,數據來源可以從inputstream,數組,集合中獲取:聚合操作可以類似SQL語句一樣的操作, 比如filter, ...
java.util.zip.ZipException: duplicate entry(重復依賴多版本的類庫)
同步SVN倉庫中的代碼,更新后,運行項目,出現如下錯誤: com.android.build.api.transform.TransformException: java.util.zip.ZipEx ...
[COGS 0065][NOIP 2002] 字串變換
65. [NOIP2002] 字串變換 ★★?? 輸入文件:string.in?? 輸出文件:string.out???簡單對比時間限制:1 s?? 內存限制:128 MB [問題描述] 已知有兩個字 ...
總結
以上是生活随笔為你收集整理的lm opencv 算法_LM算法的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 每天五分钟，玩转Docker。-Day2
- 下一篇： halcon通过点拟合圆形，鼠标选点