当前位置: 首页 > news >正文

矩阵常见分解算法及其在SLAM中的应用

文章目录

  • 常见特殊矩阵定义
  • Cholesky分解(正定Hermittian矩阵,分解结果唯一)
    • Cholesky分解应用
  • SVD分解(将singularvalues排序后分解唯一)
    • SVD 分解的应用(任意矩阵)
  • QR分解(任意矩阵,如果A可逆并且限定分解R对角线为正,则分解唯一)
    • QR分解应用

https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
在这里插入图片描述

常见特殊矩阵定义

  • 酉矩阵,正交矩阵
    • A的元素是属于复数域矩阵,如果 A A ∗ = I AA^{*} = I AA=I,那么A是属于酉矩阵(Unitary Matrix)
    • A的元素是属于实数域矩阵,如果 A A T = I AA^{T} = I AAT=I,那么A是属于正交矩阵(Orthogonal Matrix),当然也是Unitary Matrix
  • 埃尔米特矩阵,对阵矩阵
    • A的元素是属于复数域矩阵,如果 A = A ∗ A = A^{*} A=A,那么A是属于埃尔米特矩阵(Hermitian Matrix)
    • A的元素是属于实数域矩阵,如果 A = A T A = A^{T} A=AT,那么A是属于对称矩阵(Symmetric Matrix),当然也是Hermitian Matrix
  • 正定矩阵
    在这里插入图片描述

Cholesky分解(正定Hermittian矩阵,分解结果唯一)

将一个正定Hermitian矩阵分解为两个下三角矩阵的乘积
A = L L T A = LL^T A=LLT
另一种分解方式 A = L L T = L D L T A = LL^T = LDL^T A=LLT=LDLT,要求A是半正定或者半负定矩阵,条件比上面的分解方式宽松

Cholesky分解应用

  • 对于LinearGaussian 系统的状态估计问题求解 H T W − 1 H x = H T W − 1 Z H^TW^{-1}Hx = H^TW^{-1}Z HTW1Hx=HTW1Z
    求解x的时候会把 H T W − 1 H = L L ∗ H^TW^{-1}H = LL^{*} HTW1H=LL 分解完,令 L ∗ x = d L^{*}x = d Lx=d
    (1)求解 L d = H T W − 1 Z Ld = H^TW^{-1}Z Ld=HTW1Z中的d
    (2)求解 L ∗ x = d L^{*}x = d Lx=d中的x
    以上可以导出LG系统的batch递推形式,包括5个前向公式和1个后项公式。其等价于RTS算法。

SVD分解(将singularvalues排序后分解唯一)

将任意矩阵 A M × N = U M ∗ M Σ M × N V N × N A_{M \times N} = U_{M*M} \Sigma_{M \times N} V_{N \times N} AM×N=UMMΣM×NVN×N,其中 U和V都是Unitary Matrix, Σ \Sigma Σ是对角矩阵,
如果M > N,如下图。U_null_space.transpose() * A 应该是零矩阵。
反之如N> M,会有A * V_null_space为零矩阵。
MSCKF利用零空间这个性质,可以消除feature 3d位置估计误差在整个残差中的影响。

  Eigen::MatrixXd A = Eigen::MatrixXd::Random(6, 3);// Eigen::MatrixXd A = Eigen::MatrixXd::Random(3, 6);A << 1, 0, -1, -2, 1, 4, 3, 4, 5, 5, -7, 9, -1, 4, -6, 4, 7, -9;cout << "Here is the matrix A:" << endl << A << endl << endl;// do SVD decomposition, and print out the singular valuesEigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV);cout << "Its singular values are:" << endl<< svd.singularValues() << endl<< endl;cout << "Its left singular vectors are the columns of the thin U matrix:"<< endl<< svd.matrixU() << endl<< endl;cout << "Its right singular vectors are the columns of the thin V matrix:"<< endl<< svd.matrixV() << endl<< endl;const int dim = svd.singularValues().size();if (A.rows() >= A.cols()) {const Eigen::MatrixXd U_null_space =svd.matrixU().rightCols(A.rows() - dim);const Eigen::MatrixXd U_null_space_times_A = U_null_space.transpose() * A;cout << "U null space * A" << endl<< U_null_space_times_A << endl<< endl;  // should be zero} else {const Eigen::MatrixXd V_null_space =svd.matrixV().rightCols(A.cols() - dim);const Eigen::MatrixXd A_times_V_null_space = A * V_null_space;cout << "A * V null space" << endl<< A_times_V_null_space << endl<< endl;  // should be zero}cout << "U*U^T = \n"<< svd.matrixU() * svd.matrixU().transpose() << endl<< endl;  // should be identitycout << "V*V^T = \n"<< svd.matrixV() * svd.matrixV().transpose() << endl<< endl;  // should be identityreturn 0;

在这里插入图片描述

SVD 分解的应用(任意矩阵)

  • 用于PCA,SVD的奇异值的平方等于特征值,即 σ i = λ i \sigma_i = \sqrt{\lambda_i} σi=λi 。比如一堆点云,可以利用PCA性质提取line point和plane point。
  • 求矩阵零空间。MSCKF中的应用。矩阵列向量的零空间维度为:dim(null space) = cols(A) - rank(A)
    • 行向量子空间的维度 = rank(A) = m - (A的零空间的维度)= m - dim[Null(A)]
    • 列向量子空间的维度 = rank(A) = n - (A的零空间的维度)=n - dim[Null(A)]

QR分解(任意矩阵,如果A可逆并且限定分解R对角线为正,则分解唯一)

一个m*n矩阵A(m>=n),可以分解为一个Unitary Matrix Q和一个上三角矩阵R
A m × n = Q m × m R m × n A_{m \times n} = Q_{m \times m} R_{m \times n} Am×n=Qm×mRm×n
其中R的后m-n行全部为零,可以写作
A = Q R = Q [ R 1 n × n 0 ] = [ Q 1 m × n Q 2 m × ( m − n ) ] [ R 1 n × n 0 ( m − n ) × n ] = Q 1 m × n R 1 n × n A=Q R=Q\left[\begin{array}{c} R_{1_{n \times n}} \\ 0 \end{array}\right]=\left[\begin{array}{ll} Q_{1_{m \times n}} & Q_{2_{m \times (m-n)}} \end{array}\right]\left[\begin{array}{c} R_{1_{n \times n}} \\ 0_{(m-n) \times n} \end{array}\right]=Q_{1_{m \times n}} R_{1_{n \times n}} A=QR=Q[R1n×n0]=[Q1m×nQ2m×(mn)][R1n×n0(mn)×n]=Q1m×nR1n×n
类似的可以有QL,RQ,LQ分解,其中L是下三角矩阵。

  Eigen::MatrixXd A = Eigen::MatrixXd::Random(6, 3);A << 1, 0, -1, -2, 1, 4, 3, 4, 5, 5, -7, 9, -1, 4, -6, 4, 7, -9;Eigen::VectorXd b = Eigen::VectorXd::Random(6);Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(A);MatrixXd householderQ = qr.householderQ();MatrixXd matrixQ = qr.matrixQ();  // 和householderQ一样MatrixXd matrixQR_triangular_upper =qr.matrixQR().triangularView<Eigen::Upper>();MatrixXd matrixR = qr.matrixR();cout << "The rank of A is " << qr.rank() << endl << endl;cout << "householderQ matrix is:\n" << householderQ << endl << endl;cout << "matrixQ matrix is:\n" << matrixQ << endl << endl;cout << "matrixQR matrix is:\n" << qr.matrixQR() << endl << endl;cout << "matrixQR.triangularView<Eigen::Upper>() matrix is:\n"<< matrixQR_triangular_upper << endl<< endl;cout << "matrixR matrix is:\n" << matrixR << endl << endl;// Do QR decomposition to sove Ax = bEigen::VectorXd x = qr.solve(b);cout << "The solution is:\n" << x << endl << endl;

[图片]

QR分解应用

  • 求矩阵的零空间,类似SVD
  • 用来求解线性最小二乘问题 A m × n x n × 1 = b m × 1 A_{m\times n}x_{n \times 1} = b_{m \times 1} Am×nxn×1=bm×1, m >> n
    通常的解法 x = ( A T A ) m × m − 1 b x = (A^TA)_{m \times m}^{-1}b x=(ATA)m×m1b 直接求逆的话,维度很大,耗时
    将A进行QR分解,化简得到 R 1 x = Q 1 T b R_1 x = Q_1^Tb R1x=Q1Tb 并且R1是上三角矩阵,无需求逆,直接方程最后一行开始求解,便可以快速得到x
  • Kalman filter中观测量的维度m过大导致m>>状态维度n, H m × n H_{m \times n} Hm×n矩阵计算 ( H P H T + R ) − 1 (HPH^T + R)^{-1} (HPHT+R)1非常耗时
 const MatrixXX innovation_covariance = H * P * HT + V;const MatrixNX K = P * HT * innovation_covariance.inverse();

观测误差方程 r ( x ) = H X ~ + n 0 r(x) = H \tilde{X} + n_0 r(x)=HX~+n0
对H进行QR分解,两边乘以 Q T Q^T QT并带入上式得到 [ Q 1 T r o Q 2 T r o ] = [ R 1 0 ] X ~ + [ Q 1 T n o Q 2 T n o ] \left[\begin{array}{l} \mathbf{Q}_1^T \mathbf{r}_o \\ \mathbf{Q}_2^T \mathbf{r}_o \end{array}\right]=\left[\begin{array}{c} \mathbf{R}_1 \\ \mathbf{0} \end{array}\right] \tilde{\mathbf{X}}+\left[\begin{array}{l} \mathbf{Q}_1^T \mathbf{n}_o \\ \mathbf{Q}_2^T \mathbf{n}_o \end{array}\right] [Q1TroQ2Tro]=[R10]X~+[Q1TnoQ2Tno]
于是观测量维度为m的问题转化为了观测量维度为n的
Q 1 T r 0 = R 1 X ~ + Q 1 T n 0 Q_1^Tr_0 = R1 \tilde{X} + Q_1^T n_0 Q1Tr0=R1X~+Q1Tn0
将原问题中的H,残差r,观测噪声n都进行了变换,代码实现如下

  MatrixXd H_thin;VectorXd r_thin;if (H.rows() > H.cols()) {// Convert H to a sparse matrix.SparseMatrix<double> H_sparse = H.sparseView();// Perform QR decompostion on H_sparse.SPQR<SparseMatrix<double>> spqr_helper;spqr_helper.setSPQROrdering(SPQR_ORDERING_NATURAL);spqr_helper.compute(H_sparse);MatrixXd H_temp;VectorXd r_temp;(spqr_helper.matrixQ().transpose() * H).evalTo(H_temp);(spqr_helper.matrixQ().transpose() * r).evalTo(r_temp);H_thin = H_temp.topRows(21 + state_server.cam_states.size() * 6);r_thin = r_temp.head(21 + state_server.cam_states.size() * 6);// HouseholderQR<MatrixXd> qr_helper(H);// MatrixXd Q = qr_helper.householderQ();// MatrixXd Q1 = Q.leftCols(21+state_server.cam_states.size()*6);// H_thin = Q1.transpose() * H;// r_thin = Q1.transpose() * r;} else {H_thin = H;r_thin = r;}// Compute the Kalman gain.const MatrixXd &P = state_server.state_cov;MatrixXd S = H_thin * P * H_thin.transpose() +Feature::observation_noise *MatrixXd::Identity(H_thin.rows(), H_thin.rows());// MatrixXd K_transpose = S.fullPivHouseholderQr().solve(H_thin*P);MatrixXd K_transpose = S.ldlt().solve(H_thin * P);MatrixXd K = K_transpose.transpose();

相关文章:

矩阵常见分解算法及其在SLAM中的应用

文章目录 常见特殊矩阵定义Cholesky分解&#xff08;正定Hermittian矩阵&#xff0c;分解结果唯一&#xff09;Cholesky分解应用 SVD分解&#xff08;将singularvalues排序后分解唯一&#xff09;SVD 分解的应用&#xff08;任意矩阵&#xff09; QR分解&#xff08;任意矩阵&a…...

【排序】快速排序详解

✨✨欢迎大家来到Celia的博客✨✨ &#x1f389;&#x1f389;创作不易&#xff0c;请点赞关注&#xff0c;多多支持哦&#x1f389;&#x1f389; 所属专栏&#xff1a;排序 个人主页&#xff1a;Celias blog~ 一、快速排序的思想 快速排序的核心思想是&#xff1a; 选定一个…...

贪心算法总结(2)

一、买卖股票的最佳时机 . - 力扣&#xff08;LeetCode&#xff09; class Solution { public:int maxProfit(vector<int>& prices) {int miniINT_MAX;int ret0;for(int&price:prices){//遍历的时候&#xff0c;我们随时去更新最小的值&#xff0c;然后让每一位…...

弘景光电:技术实力与创新驱动并进

在光学镜头及摄像模组产品领域&#xff0c;广东弘景光电科技股份有限公司&#xff08;以下简称“弘景光电”&#xff09;无疑是一颗耀眼的明星。自成立以来&#xff0c;弘景光电凭借其强大的研发实力、卓越的产品性能、精密的制造工艺以及严格的质量管理体系&#xff0c;在光学…...

2024年7月23日~2024年7月29日周报

目录 一、前言 二、完成情况 2.1 一种具有边缘增强特点的医学图像分割网络 2.2 融合边缘增强注意力机制和 U-Net 网络的医学图像分割 2.3 遇到的困难 三、下周计划 一、前言 上周参加了一些师兄师姐的论文讨论会议&#xff0c;并完成了初稿。 本周继续修改论文&#xff0…...

M3U8流视频数据爬虫

M3U8流视频数据爬虫 HLS技术介绍 现在大部分视频客户端都采用HTTP Live Streaming&#xff08;HLS&#xff0c;Apple为了提高流播效率开发的技术&#xff09;&#xff0c;而不是直接播放MP4等视频文件。HLS技术的特点是将流媒体切分为若干【TS片段】&#xff08;比如几秒一段…...

保护您的数字财富:模块化沙箱在源代码防泄露中的突破

在数字化浪潮中&#xff0c;企业面临着前所未有的数据安全挑战。源代码、商业机密、客户数据……这些宝贵的数字资产一旦泄露&#xff0c;后果不堪设想。SDC沙盒防泄密系统&#xff0c;以其卓越的技术实力和创新的解决方案&#xff0c;为企业提供了一个坚不可摧的安全屏障。 核…...

FFmpeg源码:avio_r8、avio_rl16、avio_rl24、avio_rl32、avio_rl64函数分析

一、引言 AVIOContext是FFmpeg&#xff08;本文演示用的FFmpeg源码版本为5.0.3&#xff09;中的字节流上下文结构体&#xff0c;用来管理输入输出数据。打开一个媒体文件的时候&#xff0c;需要先把数据从硬盘读到缓冲区&#xff0c;然后会用到AVIOContext中的如下成员&#x…...

如何使用 API 查看极狐GitLab 镜像仓库中的镜像?

GitLab 是一个全球知名的一体化 DevOps 平台&#xff0c;很多人都通过私有化部署 GitLab 来进行源代码托管。极狐GitLab &#xff1a;https://gitlab.cn/install?channelcontent&utm_sourcecsdn 是 GitLab 在中国的发行版&#xff0c;专门为中国程序员服务。可以一键式部署…...

软件-vscode-plantUML-IDEA

文章目录 vscode基础命令 实操1. vscode实现springboot项目搭建 &#xff08;包括spring data jpa和sqlLite连接&#xff09; PlantUMLIDEA下载及安装Eval Reset插件配置修改IDEA创建项目的默认目录IDEA配置gitIDEA翻译插件translationIDEA断点调试IDEA全局搜索快捷键不能使用代…...

ES6语法详解,面试必会,通俗易懂版

目录 Set的基本使用WeakSet 使用Set 和 WeakSet 区别内存泄漏示例&#xff1a;使用普通 Set 保存 DOM 节点如何避免这个内存泄漏MapWeakMap 的使用 Set的基本使用 在ES6之前&#xff0c;我们存储数据的结构主要有两种&#xff1a;数组、对象。 在ES6中新增了另外两种数据结构&a…...

CTFshow--Web--代码审计

目录 web301 web302 web303 web304 web305 web306 web307 web308 web309 web310 web301 开始一个登录框, 下意识sql尝试一下 发现 1 的时候会到一个 checklogin.php 的路径下, 但啥也没有 好吧, 这是要审计代码的 ,下载好源码, 开始审计 看了一下源码 , 应该就是sql…...

Java语言程序设计——篇十(1)

&#x1f33f;&#x1f33f;&#x1f33f;跟随博主脚步&#xff0c;从这里开始→博主主页&#x1f33f;&#x1f33f;&#x1f33f; 接口介绍 接口概述接口定义接口的实现实战演练 &#x1f445;接口的继承实战演练实战演练 接口的类型常量实战演练 静态方法默认方法解决默认方…...

Qt对比MFC优势

从Qt小白到现在使用了有四年的时间&#xff0c;之前也搞过MFC,WinForm,基本上都是桌面的框架&#xff0c; 从难易程度看MFC>QT>WinForm; 运行的效率上来看MFC>QT>WinForm; 开发效率上WinForm>QT>MFC; 跨平台Qt首选&#xff1b; 界面的美观难易程度Qt>…...

RuntimeError: No CUDA GPUs are available

RuntimeError: No CUDA GPUs are available 目录 RuntimeError: No CUDA GPUs are available 【常见模块错误】 【解决方案】 解决步骤如下&#xff1a; 欢迎来到英杰社区https://bbs.csdn.net/topics/617804998 欢迎来到我的主页&#xff0c;我是博主英杰&#xff0c;211科…...

URL参数中携带中文?分享 1 段优质 JS 代码片段!

本内容首发于工粽号&#xff1a;程序员大澈&#xff0c;每日分享一段优质代码片段&#xff0c;欢迎关注和投稿&#xff01; 大家好&#xff0c;我是大澈&#xff01; 本文约 800 字&#xff0c;整篇阅读约需 1 分钟。 今天分享一段优质 JS 代码片段&#xff0c;在发送 ajax 请…...

sass的使用

一、变量 //声明一个变量 $highlight-color: #F90; .selected {border: 1px solid $highlight-color; }//编译后 .selected {border: 1px solid #F90; }二、导入 import "xxx.scss"三、混合器简单定义 通过mixin定义&#xff0c;通过include调用 // mixin.scss /…...

【足球走地软件】走地数据分析预测【大模型篇】走地预测软件实战分享

了解什么是走地数据&#xff1f; 走地数据分析&#xff0c;在足球赛事的上下文中&#xff0c;是一种针对正在进行中的比赛进行实时数据分析的方法。这种方法主要用于预测比赛中的某些结果或趋势&#xff0c;如总进球数、比分变化、球队表现等。 在足球走地数据分析中&#xf…...

现在有什么赛道可以干到退休?

最近&#xff0c;一则“90后无论男女都得65岁以后退休”的消息在多个网络平台流传&#xff0c;也不知道是真是假&#xff0c;好巧不巧今天刷热点的时候又看到一条这样的热点&#xff1a;现在有什么赛道可以干到退休&#xff1f; 点进去看了几条热评&#xff0c;第一条热评说的…...

c程序杂谈系列(职责链模式与if_else)

从处理器的角度来说&#xff0c;条件分支会导致指令流水线的中断&#xff0c;所以控制语句需要严格保存状态&#xff0c;因为处理器是很难直接进行逻辑判断的&#xff0c;有可能它会执行一段时间&#xff0c;发现出错后再返回&#xff0c;也有可能通过延时等手段完成控制流的正…...

进程地址空间(比特课总结)

一、进程地址空间 1. 环境变量 1 &#xff09;⽤户级环境变量与系统级环境变量 全局属性&#xff1a;环境变量具有全局属性&#xff0c;会被⼦进程继承。例如当bash启动⼦进程时&#xff0c;环 境变量会⾃动传递给⼦进程。 本地变量限制&#xff1a;本地变量只在当前进程(ba…...

五年级数学知识边界总结思考-下册

目录 一、背景二、过程1.观察物体小学五年级下册“观察物体”知识点详解&#xff1a;由来、作用与意义**一、知识点核心内容****二、知识点的由来&#xff1a;从生活实践到数学抽象****三、知识的作用&#xff1a;解决实际问题的工具****四、学习的意义&#xff1a;培养核心素养…...

【项目实战】通过多模态+LangGraph实现PPT生成助手

PPT自动生成系统 基于LangGraph的PPT自动生成系统&#xff0c;可以将Markdown文档自动转换为PPT演示文稿。 功能特点 Markdown解析&#xff1a;自动解析Markdown文档结构PPT模板分析&#xff1a;分析PPT模板的布局和风格智能布局决策&#xff1a;匹配内容与合适的PPT布局自动…...

解决本地部署 SmolVLM2 大语言模型运行 flash-attn 报错

出现的问题 安装 flash-attn 会一直卡在 build 那一步或者运行报错 解决办法 是因为你安装的 flash-attn 版本没有对应上&#xff0c;所以报错&#xff0c;到 https://github.com/Dao-AILab/flash-attention/releases 下载对应版本&#xff0c;cu、torch、cp 的版本一定要对…...

BCS 2025|百度副总裁陈洋:智能体在安全领域的应用实践

6月5日&#xff0c;2025全球数字经济大会数字安全主论坛暨北京网络安全大会在国家会议中心隆重开幕。百度副总裁陈洋受邀出席&#xff0c;并作《智能体在安全领域的应用实践》主题演讲&#xff0c;分享了在智能体在安全领域的突破性实践。他指出&#xff0c;百度通过将安全能力…...

网络编程(UDP编程)

思维导图 UDP基础编程&#xff08;单播&#xff09; 1.流程图 服务器&#xff1a;短信的接收方 创建套接字 (socket)-----------------------------------------》有手机指定网络信息-----------------------------------------------》有号码绑定套接字 (bind)--------------…...

视频行为标注工具BehaviLabel(源码+使用介绍+Windows.Exe版本)

前言&#xff1a; 最近在做行为检测相关的模型&#xff0c;用的是时空图卷积网络&#xff08;STGCN&#xff09;&#xff0c;但原有kinetic-400数据集数据质量较低&#xff0c;需要进行细粒度的标注&#xff0c;同时粗略搜了下已有开源工具基本都集中于图像分割这块&#xff0c…...

【从零学习JVM|第三篇】类的生命周期(高频面试题)

前言&#xff1a; 在Java编程中&#xff0c;类的生命周期是指类从被加载到内存中开始&#xff0c;到被卸载出内存为止的整个过程。了解类的生命周期对于理解Java程序的运行机制以及性能优化非常重要。本文会深入探寻类的生命周期&#xff0c;让读者对此有深刻印象。 目录 ​…...

08. C#入门系列【类的基本概念】:开启编程世界的奇妙冒险

C#入门系列【类的基本概念】&#xff1a;开启编程世界的奇妙冒险 嘿&#xff0c;各位编程小白探险家&#xff01;欢迎来到 C# 的奇幻大陆&#xff01;今天咱们要深入探索这片大陆上至关重要的 “建筑”—— 类&#xff01;别害怕&#xff0c;跟着我&#xff0c;保准让你轻松搞…...

【Linux】Linux 系统默认的目录及作用说明

博主介绍&#xff1a;✌全网粉丝23W&#xff0c;CSDN博客专家、Java领域优质创作者&#xff0c;掘金/华为云/阿里云/InfoQ等平台优质作者、专注于Java技术领域✌ 技术范围&#xff1a;SpringBoot、SpringCloud、Vue、SSM、HTML、Nodejs、Python、MySQL、PostgreSQL、大数据、物…...