- Java锁的逻辑(结合对象头和ObjectMonitor)
- 还在用饼状图?来瞧瞧这些炫酷的百分比可视化新图形(附代码实现)⛵
- 自动注册实体类到EntityFrameworkCore上下文,并适配ABP及ABPVNext
- 基于Sklearn机器学习代码实战
\(\quad\) 找到一个 n 维的变量 \(\mathbf{x}^{*} \in \mathbb{R}^{n}\) , 使得损失函数 \(F(\mathbf{x})\) 取局部最小值:
\(\quad\) 其中 \(f_{i}\) 是残差函数, 比如测量值和预测值之间的差, 且有 \(m \geq n\) 。 部最小值指对任意 \(\left\|\mathbf{x}-\mathbf{x}^{*}\right\|<\delta\) 有 \(F\left(\mathbf{x}^{*}\right) \leq F(\mathbf{x})\) \(\quad\) 损失函数泰勒展开,假设损失函数 \(F(\mathbf{x})\) 是可导并且平滑的, 因此, 二阶泰勒展开
这里要着重注意一下,这里的 \(\mathbf{J}\) 和 \(\mathbf{H}\) 都是每一个残差项的雅可比堆叠(计算)而来,实际上对于初学者来说并不直观,后面我们会以一个曲线拟合和 \(BA\) 问题来详细分析一下如何通过连加来推算到 \(\mathbf{J}\) 和 \(\mathbf{H}\) 。
\(\quad\) 其中 \(\mathbf{J}\) 和 \(\mathbf{H}\) 分别为损失函数 \(F\) 对变量 \(x\) 的一阶导和二阶导矩阵,也就是我们通常所说的雅可比矩阵和海塞矩阵。( 这里的 \(\mathbf{x}\) 包含了所有待优化的变量,在视觉SLAM问题中,一般是相机的 Pose 和已经三角化的点的坐标或者逆深度,且由于相机一般能观测到的3D点的个数是有限的,因此其雅可比矩阵也就是稀疏的,只有两个地方的雅可比求导是不为0的,参考14讲P247,那么 \(J_{i,j}^TJ_{i,j}\) ,则只有四个地方是不为0的).
\(\quad\) 忽略泰勒展开的高阶项,损失函数变成了二次函数,可以轻易得到如下性质:
迭代法初衷 :
找到一个下降方向使得损失函数随着 \(x\) 的迭代逐渐减少直到 \(x^*\) .
分两个步骤;第一,找到下降方向 单位向量 \(d\) ,第二,确定下降的 步长 \(a\) .
假设 \(a\) 足够的小,又因为 \(d\) 为单位向量,因此可以将 \(ad\) 看作是一个微小的变化量 \(\triangle{x}\) ,我们可以对损失函数进行一阶泰勒展开:
只需要寻找下降方向,满足:
通过 line search 的方法找到下降的步长: \(a^*=argmin_{a>0} [F(x+ad)]\) 。
适用于迭代的开始阶段 。
\(\quad\) 从下降方向的条件(单位向量)可以知道: \(\mathbf{Jd=||J||}cos\theta\) ,其中 \(\theta\) 表示的是下降方向和梯度方向的夹角. 当 \(\theta = \pi\) 有
这里为什么能写成向量的内积运算,笔者在刚开始看起来还认为是两个矩阵相乘法,却直接写成了内积运算,仔细思索发现 \(d\) 其实上是一个和 \(x\) 相同维度的单位向量,其纬度为 \(n\times 1\) ,而雅可比矩阵 。
\(\quad\) 即 梯度的负方向 为最速下降方向。缺点:最优值附近震荡,收敛慢.
在局部最优点 \(x^∗\) 附近,如果 \(x + ∆x\) 是最优解,则损失函数对 \(∆x\) 的导数等于 \(0\) ,对公式 (2) 取一阶导有
得到: \(∆x = -\mathbf{H^{-1}J^T}\) 。缺点:二阶导矩阵计算复杂.
这里我们其实既可以看作是多个残差的分量相加后组成的 \(H\) ,也可以看作是每个残差单独的 \(H\) .
将损失函数的二阶泰勒展开记作:
求以下函数的最小化
其中, \(μ ≥ 0\) 为阻尼因子, $ \frac{1}{2}\mu \Delta x^T\Delta x $是惩罚项。对新的损失函数求一阶导,并令其等于 \(0\) 有
残差函数 \(f(x)\) 为非线性函数,对其进行一阶泰勒近似有:
带入损失函数:
这里我们假设暂时只关注其中的一项(其实也可以是所有损失项的叠加,最终形式是一样的)。在 \(x\) 处进行的泰勒展开,则认为 \(x\) 是已知的,现在的损失函数是一个关于 \(\Delta x\) 的函数,其最小值,则令关于 \(\Delta x\) 的导数为 \(0\) 即可。可以得到:
上面这个形式就是我们在论文或者各种SLAM问题中经常见到的形式了, \(\mathbf{H \Delta x =b}\) ,也叫做 normal equation 。
现在我们通过非线性最小二乘来进行线性拟合实验,将理论应用于实际中去。假设曲线方程为:
其中设 \(a=1,b=2,c=3\) .
现在加入噪声项生成带有高斯分布的噪声数据,当然不是高斯分布的数据也是可以的,但是在自己实验的时候尽量不要出现外点数据,因为我们并没有处理外点数据的策略。则生成数据的方程为:
其中 \(w\) 为符合高斯分布的噪声数据.
通过如下程序生成观测数据:
double ar = 1.0, br = 2.0, cr = 1.0; // 真实参数值
int N = 100; // 数据点
double w_sigma = 1.0; // 噪声Sigma值
vector<double> x_data, y_data; // 数据
for (int i = 0; i < N; i++) {
double x = i / 100.0;
x_data.push_back(x);
y_data.push_back(exp(ar * x * x + br * x + cr) + rng.gaussian(w_sigma * w_sigma));
}
接下来我们关心雅可比如何计算,误差项 \(f_i(a,b,c)\) 可以写成如下形式:
我们知道这两项相减是绝对不可能相等的,因为在生成数据的时候加入了高斯噪声。我们这里有 \(N\) 个观测,即 \(i\in (1-100)\) ,我们将其写成连加的形式 。
该式中右边就是残差项的具体形式,我们对其进行平方,防止负的残差和正的残差抵消的情况,前面我们已经说过可以将残差项通过一阶泰勒展开进行近似,然后进行平方得到每一个残差项的具体形式:
\(f(x+\Delta x)\approx f(x)+J(x)\Delta x\) 。
此时由于某时刻的观测已知,因此误差项是一个关于 \(\Delta x\) 的二次函数,求该项的最小值只要让关于 \(\Delta x\) 的导数为 \(0\) 即可。求导后可得:
这里我们简单的记:
\[\boldsymbol{J}(\boldsymbol{x})^{T} f(\boldsymbol{x}) = \mathbf{\eta}\\ \boldsymbol{J}(\boldsymbol{x})^{T} \boldsymbol{J}(\boldsymbol{x}) \Delta \boldsymbol{x}=\mathbf{H\Delta x} \]
即我们常见的形式:
读者要注意到这里的 \(b\) 其实就是上面的 \(-\eta\) 。
这里我们假设残差项记为 \(\mathbf{e_i}\) 一共有 \(N\) 个观测,则有 \(N\) 个残差项.
整个 \(F(X)\) 此时是关于待优化变量的函数,每一项分别用各自的一阶泰勒展开近似,注意这里的每一项由于观测的不同,每一项都是一个不同的函数表达式,但是优化变量都是一样的。得到如下结果:
\[\begin{aligned} \frac{1}{2}\|f(\boldsymbol{x})+\boldsymbol{J}(\boldsymbol{x}) \Delta \boldsymbol{x}\|^{2} &=\Omega(\Delta x) \end{aligned} \]
这里的 \(\Delta x\) 是我们在使用基于迭代下降的方法中所选中的步长和方向,如果 \(F(X)\) 在 \(\Delta x\) 为某个值时取得极小值,则 \(\Delta x\) 无论是在任何一个方向加或者减函数值都会上升,此时这个点则为极小值点,这里的叙述不太数学化,但是大家联想一下极小值的定义,应该是可以理解的,当达到该条件后,那么该点关于 \(\Delta x\) 的导数一定为 \(0\) 。所以对此时的 \(F(X)\) 求导并让其等于 \(0\) 得到:
再将该式子变形,将关于 \(\Delta x\) 的项都移动到左边,没有关于 \(\Delta x\) 的移动到右边:
其实也就是
写成连加的形式:
这里我们就通过每一项的一个具体形式来推倒出最后的 H 和 b 是怎么来的了。也就是我们经常在程序中见到的 += 操作的原理:
H += J * J.transpose();
b += -J * error;
我们再次回到曲线拟合的题目中去,待优化的变量就三个 \(a,b,c\) 则每一个残差项都含有这三个参数,我们称其雅可比为稠密的(虽然只有三个参数,视觉BA问题中由于相机观测的特殊性,其雅可比矩阵是稀疏的),对每一个残差向分别求雅可比,然后求和得到最终的 \(H\) 和 \(b\) ,然后求解一次 \(\Delta x\) ,Step 一次,根据判断条件选择是否继续进行迭代。每一个残差项对于 \(\Delta x\) 的雅可比为 。
得到了雅可比,那么剩下的就是迭代求解即可,完整代码如下,来自14讲配套代码:
#include <iostream>
#include <chrono>
#include <opencv2/opencv.hpp>
#include <Eigen/Core>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main(int argc, char **argv) {
double ar = 1.0, br = 2.0, cr = 1.0; // 真实参数值
double ae = 2.0, be = -1.0, ce = 5.0; // 估计参数值
int N = 100; // 数据点
double w_sigma = 1.0; // 噪声Sigma值
double inv_sigma = 1.0 / w_sigma;
cv::RNG rng; // OpenCV随机数产生器
vector<double> x_data, y_data; // 数据
for (int i = 0; i < N; i++) {
double x = i / 100.0;
x_data.push_back(x);
y_data.push_back(exp(ar * x * x + br * x + cr) + rng.gaussian(w_sigma * w_sigma));
}
// 开始Gauss-Newton迭代
int iterations = 100; // 迭代次数
double cost = 0, lastCost = 0; // 本次迭代的cost和上一次迭代的cost
chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
for (int iter = 0; iter < iterations; iter++) {
Matrix3d H = Matrix3d::Zero(); // Hessian = J^T W^{-1} J in Gauss-Newton
Vector3d b = Vector3d::Zero(); // bias
cost = 0;
for (int i = 0; i < N; i++) {
double xi = x_data[i], yi = y_data[i]; // 第i个数据点
double error = yi - exp(ae * xi * xi + be * xi + ce);
Vector3d J; // 雅可比矩阵
J[0] = -xi * xi * exp(ae * xi * xi + be * xi + ce); // de/da
J[1] = -xi * exp(ae * xi * xi + be * xi + ce); // de/db
J[2] = -exp(ae * xi * xi + be * xi + ce); // de/dc
H += J * J.transpose();
b += -J * error;
cost += error * error;
}
// 求解线性方程 Hx=b
Vector3d dx = H.ldlt().solve(b);
if (isnan(dx[0])) {
cout << "result is nan!" << endl;
break;
}
if (iter > 0 && cost >= lastCost) {
cout << "cost: " << cost << ">= last cost: " << lastCost << ", break." << endl;
break;
}
ae += dx[0];
be += dx[1];
ce += dx[2];
lastCost = cost;
cout << "total cost: " << cost << ", \t\tupdate: " << dx.transpose() <<
"\t\testimated params: " << ae << "," << be << "," << ce << endl;
}
chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
chrono::duration<double> time_used = chrono::duration_cast<chrono::duration<double>>(t2 - t1);
cout << "solve time cost = " << time_used.count() << " seconds. " << endl;
cout << "estimated abc = " << ae << ", " << be << ", " << ce << endl;
return 0;
}
total cost: 3.19575e+06, update: 0.0455771 0.078164 -0.985329 estimated params: 2.04558,-0.921836,4.01467
total cost: 376785, update: 0.065762 0.224972 -0.962521 estimated params: 2.11134,-0.696864,3.05215
total cost: 35673.6, update: -0.0670241 0.617616 -0.907497 estimated params: 2.04432,-0.0792484,2.14465
total cost: 2195.01, update: -0.522767 1.19192 -0.756452 estimated params: 1.52155,1.11267,1.3882
total cost: 174.853, update: -0.537502 0.909933 -0.386395 estimated params: 0.984045,2.0226,1.00181
total cost: 102.78, update: -0.0919666 0.147331 -0.0573675 estimated params: 0.892079,2.16994,0.944438
total cost: 101.937, update: -0.00117081 0.00196749 -0.00081055 estimated params: 0.890908,2.1719,0.943628
total cost: 101.937, update: 3.4312e-06 -4.28555e-06 1.08348e-06 estimated params: 0.890912,2.1719,0.943629
total cost: 101.937, update: -2.01204e-08 2.68928e-08 -7.86602e-09 estimated params: 0.890912,2.1719,0.943629
cost: 101.937>= last cost: 101.937, break.
solve time cost = 0.00440302 seconds.
estimated abc = 0.890912, 2.1719, 0.943629
\(a\) | \(b\) | \(c\) | |
---|---|---|---|
Estimate | \(0.890912\) | \(2.1719\) | \(0.943629\) |
Real | \(1\) | \(2\) | \(1\) |
下一节我们来讨论一下视觉SLAM中的非线性优化问题的具体形式,以及其 \(H\) 和 \(b\) 的由来和构建方法.
最后此篇关于非线性优化问题基本形式概述的文章就讲到这里了,如果你想了解更多关于非线性优化问题基本形式概述的内容请搜索CFSDN的文章或继续浏览相关文章,希望大家以后支持我的博客! 。
关闭。这个问题是off-topic .它目前不接受答案。 想要改进这个问题? Update the question所以它是on-topic用于堆栈溢出。 关闭 12 年前。 Improve thi
我有一个动态网格,其中的数据功能需要正常工作,这样我才能逐步复制网格中的数据。假设在第 5 行中,我输入 10,则从第 6 行开始的后续行应从 11 开始读取,依此类推。 如果我转到空白的第一行并输入
我有一个关于我的按钮消失的问题 我已经把一个图像作为我的按钮 用这个函数动画 function example_animate(px) { $('#cont
我有一个具有 Facebook 连接和经典用户名/密码登录的网站。目前,如果用户单击 facebook_connect 按钮,系统即可运行。但是,我想将现有帐户链接到 facebook,因为用户可以选
我有一个正在为 iOS 开发的应用程序,该应用程序执行以下操作 加载和设置注释并启动核心定位和缩放到位置。 map 上有很多注释,从数据加载不会花很长时间,但将它们实际渲染到 map 上需要一段时间。
我被推荐使用 Heroku for Ruby on Rails 托管,到目前为止,我认为我真的会喜欢它。只是想知道是否有人可以帮助我找出问题所在。 我按照那里的说明在该网站上创建应用程序,创建并提交
我看过很多关于 SSL 错误的帖子和信息,我自己也偶然发现了一个。 我正在尝试使用 GlobalSign CA BE 证书通过 Android WebView 访问网页,但出现了不可信错误。 对于大多
我想开始使用 OpenGL 3+ 和 4,但我在使用 Glew 时遇到了问题。我试图将 glew32.lib 包含在附加依赖项中,并且我已将库和 .dll 移动到主文件夹中,因此不应该有任何路径问题。
我已经盯着这两个下载页面的源代码看了一段时间,但我似乎找不到问题。 我有两个下载页面,一个 javascript 可以工作,一个没有。 工作:http://justupload.it/v/lfd7不是
我一直在使用 jQuery,只是尝试在单击链接时替换文本字段以及隐藏/显示内容项。它似乎在 IE 中工作得很好,但我似乎无法让它在 FF 中工作。 我的 jQuery: $(function() {
我正在尝试为 NDK 编译套接字库,但出现以下两个错误: error: 'close' was not declared in this scope 和 error: 'min' is not a m
我正在使用 Selenium 浏览器自动化框架测试网站。在测试过程中,我切换到特定的框架,我们将其称为“frame_1”。后来,我在 Select 类中使用了 deselectAll() 方法。不久之
我正在尝试通过 Python 创建到 Heroku PostgreSQL 数据库的连接。我将 Windows10 与 Python 3.6.8 和 PostgreSQL 9.6 一起使用。 我从“ht
我有一个包含 2 列的数据框,我想根据两列之间的比较创建第三列。 所以逻辑是:第 1 列 val = 3,第 2 列 val = 4,因此新列值什么都没有 第 1 列 val = 3,第 2 列 va
我想知道如何调试 iphone 5 中的 css 问题。 我尝试使用 firelite 插件。但是从纵向旋转到横向时,火石占据了整个屏幕。 有没有其他方法可以调试 iphone 5 中的 css 问题
所以我有点难以理解为什么这不起作用。我正在尝试替换我正在处理的示例站点上的类别复选框。我试图让它做以下事情:未选中时以一种方式出现,悬停时以另一种方式出现(选中或未选中)选中时以第三种方式出现(而不是
Javascript CSS 问题: 我正在使用一个文本框来写入一个 div。我使用以下 javascript 获取文本框来执行此操作: function process_input(){
你好,我很难理解 P、NP 和多项式时间缩减的主题。我试过在网上搜索它并问过我的一些 friend ,但我没有得到任何好的答案。 我想问一个关于这个话题的一般性问题: 设 A,B 为 P 中的语言(或
你好,我一直在研究 https://leetcode.com/problems/2-keys-keyboard/并想到了这个动态规划问题。 您从空白页上的“A”开始,完成后得到一个数字 n,页面上应该
我正在使用 Cocoapods 和 KIF 在 Xcode 服务器上运行持续集成。我已经成功地为一个项目设置了它来报告每次提交。我现在正在使用第二个项目并收到错误: Bot Issue: warnin
我是一名优秀的程序员,十分优秀!