lm opencv 算法_LM算法
最小二乘法的概念
最小二乘法要關(guān)心的是對應(yīng)的cost function是線性還是非線性函數(shù)，不同的方法計(jì)算效率如何，要不要求逆，矩陣的維數(shù)
一般都是過約束，方程式的數(shù)目多于未知的參數(shù)數(shù)目。
最小二乘法的目標(biāo)：求誤差的最小平方和，根據(jù)cost function的對應(yīng)有兩種：線性和非線性(取決于對應(yīng)的殘差(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\) 的維數(shù)是 1*1(y的列數(shù)目).所以 \(b^TAx = x^TA^Tb\)
\(S(x) = b^Tb - 2b^TAx + x^TA^TAx\)
對\(x\)求導(dǎo)
\(-2A^Tb + 2 (A^TA)x = 0\)
\(x = (A^TA)^{-1}A^Tb\)
而非線性最小二乘沒有closed-form，通常用迭代法求解。在每次迭代的過程使用一個線性化的方程代替計(jì)算。
有牛頓法，牛頓高斯法，LM, 其實(shí)可以分為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\)是泰勒展開式的點(diǎn)， \(\alpha\)是 步長(一個實(shí)數(shù)) , $ \overrightarrow d$ 單位方向(一個向量)，即 \(|\overrightarrow d| = 1\)
顯然當(dāng)\(d\)的方向和\(f'(x_0)\)方向相反的時候，\(cos180 =-1\),整體取到最小值。這就是為什么取的是梯度的反方向原因。
當(dāng)然上面的步長一般也是可以求得，怎么求步長也是一門學(xué)問。
下面的都是從牛頓法引申出來的,記住牛頓法求得是穩(wěn)定點(diǎn)\(f'(x) = 0\),導(dǎo)數(shù)為0的不一定是最小值,梯度下降法求得是局部最小值,從計(jì)算上看
牛頓法 泰勒展開式到二階，
\[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)
\]
求導(dǎo)就有，\(\frac{\partial f}{\partial x_t} = g + H(x_{t+1} -x_t)\)，讓他為0，就有了牛頓公式
\[x_{t+1} = x_t - H^{-1}g
\]
要是H是正定的，上面的就是凸函數(shù)，也就一定有了最小值。可惜H不一定是正定的，這就引導(dǎo)出了下面的方法
高斯-牛頓法
是另一種經(jīng)常用于求解非線性最小二乘的迭代法(一定程度上可視為標(biāo)準(zhǔn)非線性最小二乘求解方法)。
我們優(yōu)化的cost function是 $$min \sum _i r_i(x)^2$$
根據(jù)上面的牛頓法：
\[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})
\]
要是把上面公式的最后一項(xiàng)去掉，至少是半正定了，而且不用計(jì)算Hessain矩陣了，這就是牛頓高斯法
\[H_{jk} \approx 2 \sum _i J_{ij}J_{ik} \quad With \quad J_{ij} = \frac{\partial r_i}{\partial x_j}
\]
在什么情況下，去掉最后一項(xiàng)比較好，
residual(\(r_i\))小的時候
或者接近linear ,這樣一階微分是常數(shù)，二階微分就是0
上面兩種情況下，第二項(xiàng)都很小,公式如下
\[x_{t+1} = x_t - (J^TJ)^{-1}J^Tr
\]
Levenberg-Marquardt
的迭代法用于求解非線性最小二乘問題，就結(jié)合了梯度下降和高斯-牛頓法。
\[x_{t+1} = x_t - (H + \lambda I_n)^{-1}g
\]
\[x_{t+1} = x_t - (J^TJ + \lambda I_n)^{-1}J^Tr
\]
總結(jié)
所以如果把最小二乘看做是優(yōu)化問題的話，那么梯度下降是求解方法的一種，\(x=(A^TA)^{-1}A^Tb\)是求解線性最小二乘的一種，高斯-牛頓法和Levenberg-Marquardt則能用于求解非線性最小二乘。
LM算法相對于高斯牛頓算法和梯度下降的優(yōu)缺點(diǎn)
首先梯度下降法和高斯牛頓法都是最優(yōu)化方法。其區(qū)別之處在于，
梯度下降法在尋找目標(biāo)函數(shù)極小值時，是沿著反梯度方向進(jìn)行尋找的。梯度的定義就是指向標(biāo)量場增長最快的方向，在尋找極小值時，先隨便定初始點(diǎn)(x0，y0)然后進(jìn)行迭代不斷尋找直到梯度的模達(dá)到預(yù)設(shè)的要求。但是梯度下降法的缺點(diǎn)之處在于：在遠(yuǎn)離極小值的地方下降很快，而在靠近極小值的地方下降很慢，靠近的時候可能成zig-zag下降。
而高斯牛頓法是一種非線性最小二乘最優(yōu)化方法。其利用了目標(biāo)函數(shù)的泰勒展開式把非線性函數(shù)的最小二乘化問題化為每次迭代的線性函數(shù)的最小二乘化問題。高斯牛頓法的缺點(diǎn)在于：若初始點(diǎn)距離極小值點(diǎn)過遠(yuǎn)，迭代步長過大會導(dǎo)致迭代下一代的函數(shù)值不一定小于上一代的函數(shù)值。
LM算法在高斯牛頓法中加入了因子μ，當(dāng)μ大時相當(dāng)于梯度下降法，μ小時相當(dāng)于高斯牛頓法。在使用Levenberg-Marquart時，先設(shè)置一個比較小的μ值，當(dāng)發(fā)現(xiàn)目標(biāo)函數(shù)反而增大時，將μ增大使用梯度下降法快速尋找，然后再將μ減小使用牛頓法進(jìn)行尋找。
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.
% 計(jì)算函數(shù)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算法
% 初始猜測初始點(diǎn)
a0=100; b0=100;
y_init = a0*cos(b0*data_1) + b0*sin(a0*data_1);
% 數(shù)據(jù)個數(shù)
Ndata=length(obs_1);
% 參數(shù)維數(shù)
Nparams=2;
% 迭代最大次數(shù)
n_iters=60;
% LM算法的阻尼系數(shù)初值
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ù)當(dāng)前估計(jì)值，計(jì)算雅克比矩陣
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
% 根據(jù)當(dāng)前參數(shù)，得到函數(shù)值
y_est = a_est*cos(b_est*data_1) + b_est*sin(a_est*data_1);
% 計(jì)算誤差
d=obs_1-y_est;
% 計(jì)算(擬)海塞矩陣
H=J'*J;
% 若是第一次迭代，計(jì)算誤差
if it==1
e=dot(d,d);
end
end
% 根據(jù)阻尼系數(shù)lamda混合得到H矩陣
H_lm=H+(lamda*eye(Nparams,Nparams));
% 計(jì)算步長dp，并根據(jù)步長計(jì)算新的可能的\參數(shù)估計(jì)值
dp=inv(H_lm)*(J'*d(:));
%求誤差大小
g = J'*d(:);
a_lm=a_est+dp(1);
b_lm=b_est+dp(2);
% 計(jì)算新的可能估計(jì)值對應(yīng)的y和計(jì)算殘差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);
% 根據(jù)誤差，決定如何更新參數(shù)和阻尼系數(shù)
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
%顯示優(yōu)化的結(jié)果
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;//估計(jì)誤差
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)//成功則更新向量與估計(jì)誤差
{
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 };
//初始點(diǎn)
double p0[] = { 100, 100 };
LM(p0, 2, data, 32, obs, 0.1, 10);
system("pause");
}
LM算法與非線性最小二乘問題
摘錄的一篇有關(guān)求解非線性最小二乘問題的算法--LM算法的文章,當(dāng)中也加入了一些我個人在求解高精度最小二乘問題時候的一些感觸: LM算法,全稱為Levenberg-Marquard算法,它可用于解決非線 ...
Levenberg-Marquardt迭代(LM算法)-改進(jìn)Guass-Newton法
1.前言 ?????????????????????????????? a.對于工程問題,一般描述為:從一些測量值(觀測量)x 中估計(jì)參數(shù) p?即x = f(p), ??? ...
梯度下降法、牛頓法、高斯牛頓法、LM最優(yōu)化算法
1.梯度下降法 2.牛頓法 3.高斯牛頓法 4.LM算法
Levenberg-Marquardt優(yōu)化算法以及基于LM的BP-ANN
一.LM最優(yōu)化算法 ? ? 最優(yōu)化是尋找使得目標(biāo)函數(shù)有最大或最小值的的參數(shù)向量.根據(jù)求導(dǎo)數(shù)的方法,可分為2大類.(1)若f具有解析函數(shù)形式,知道x后求導(dǎo)數(shù)速度快.(2)使用數(shù)值差分來求導(dǎo)數(shù).根據(jù)使用模 ...
LM擬合算法
一.? Levenberg-Marquardt算法 (1)y=a*e.^(-b*x)形式擬合 clear all % 計(jì)算函數(shù)f的雅克比矩陣,是解析式 syms a b y x real; f=a*e ...
Levenberg-Marquardt算法基礎(chǔ)知識
Levenberg-Marquardt算法基礎(chǔ)知識 (2013-01-07 16:56:17) 轉(zhuǎn)載▼ ? 什么是最優(yōu)化?Levenberg-Marquardt算法是最優(yōu)化算法中的一種.最優(yōu)化是尋找使 ...
相機(jī)標(biāo)定：關(guān)于用Levenberg-Marquardt算法在相機(jī)標(biāo)定中應(yīng)用
LM算法在相機(jī)標(biāo)定的應(yīng)用共有三處. (1)單目標(biāo)定或雙目標(biāo)定中,在內(nèi)參固定的情況下,計(jì)算最佳外參.OpenCV中對應(yīng)的函數(shù)為findExtrinsicCameraParams2. (2)單目標(biāo)定中,在 ...
點(diǎn)云匹配和ICP算法概述
Iterative Closest Point (ICP) [1][2][3] is an algorithm employed to minimize the difference between ...
隨機(jī)推薦
Hibernate中事務(wù)聲明
Hibernate中JDBC事務(wù)聲明,在Hibernate配置文件中加入如下代碼,不做聲明Hibernate默認(rèn)就是JDBC事務(wù). 一個JDBC 不能跨越多個數(shù)據(jù)庫. Hibernate中JTA事務(wù)聲 ...
Redis基礎(chǔ)知識之————空間換時間的查詢案例
空間與時間 空間換時間是在數(shù)據(jù)庫中經(jīng)常出現(xiàn)的術(shù)語,簡單說就是把查詢需要的條件進(jìn)行索引的存儲,然后查詢時為O(1)的時間復(fù)雜度來快速獲取數(shù)據(jù),從而達(dá)到了使用空間存儲來換快速的時間響應(yīng)!對于redis這個 ...
Basic Printing Architecture
https://blogs.technet.microsoft.com/askperf/2007/06/19/basic-printing-architecture/ Printer sharing, ...
Android 文字繪制(DrawText)技術(shù)總結(jié)
這里的繪制文字不是直接調(diào)用TextView.setText(String content)去展示文字內(nèi)容.而是在View上面通過?canvas.drawText(text, x, y,textPain ...
iOS學(xué)習(xí)——如何在mac上獲取開發(fā)使用的模擬器的資源以及模擬器中每個應(yīng)用的應(yīng)用沙盒
如題,本文主要研究如何在mac上獲取開發(fā)使用的模擬器的資源以及模擬器中每個應(yīng)用的應(yīng)用沙盒.做過安卓開發(fā)的小伙伴肯定很方便就能像打開資源管理器一樣查看我們寫到手機(jī)本地或應(yīng)用中的各種資源,但是在iOS開發(fā) ...
洛谷P3835 【模板】可持久化平衡樹
題目背景 本題為題目?普通平衡樹?的可持久化加強(qiáng)版. 數(shù)據(jù)已經(jīng)經(jīng)過強(qiáng)化 感謝@Kelin 提供的一組hack數(shù)據(jù) 題目描述 您需要寫一種數(shù)據(jù)結(jié)構(gòu)(可參考題目標(biāo)題),來維護(hù)一些數(shù),其中需要提供以下操作( ...
Stream初步應(yīng)用
一.什么是stream Stream(流)是一個來自數(shù)據(jù)源的元素隊(duì)列并支持聚合操作,數(shù)據(jù)來源可以從inputstream,數(shù)組,集合中獲取:聚合操作可以類似SQL語句一樣的操作, 比如filter, ...
java.util.zip.ZipException: duplicate entry(重復(fù)依賴多版本的類庫)
同步SVN倉庫中的代碼,更新后,運(yùn)行項(xiàng)目,出現(xiàn)如下錯誤: com.android.build.api.transform.TransformException: java.util.zip.ZipEx ...
[COGS 0065][NOIP 2002] 字串變換
65. [NOIP2002] 字串變換 ★★?? 輸入文件:string.in?? 輸出文件:string.out???簡單對比時間限制:1 s?? 內(nèi)存限制:128 MB [問題描述] 已知有兩個字 ...
總結(jié)
以上是生活随笔為你收集整理的lm opencv 算法_LM算法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 每天五分钟，玩转Docker。-Day2
- 下一篇： halcon通过点拟合圆形，鼠标选点