奇异值分解（SVD）在基于ARM Cortex-M架构的嵌入式系统中的实现
摘要
奇异值分解（SVD）是嵌入式系统中若干新兴应用的核心数学工具，例如多输入多输出（MIMO）系统、数据分析、信号稀疏表示。由于SVD算法归结为求解一个计算成本高的特征值问题，已有研究提出专用硬件解决方案和并行实现以克服这一瓶颈。然而，这些方案通常需要额外的硬件资源，而这类资源在嵌入式系统中往往不可用，因此在此背景下迫切需要优化算法。本文旨在提出一种在ARM Cortex-M上高效实现的SVD算法。为此，我们进行了以下工作：（i）全面阐述最常用的SVD算法，采用统一的符号体系，提供对这些算法较为完整且深入的综述；（ii）在ARM Cortex-M4F微控制器上实现这些算法，以开发适用于无操作系统的嵌入式系统的库；（iii）通过所提出SVD算法的比较研究，找出最适合低资源裸机嵌入式系统的最佳实现方案；（iv）展示在惯性测量单元（IMU）的卡尔曼滤波中的实际应用，作为SVD如何提升现有算法精度以及其在低资源系统中实用性的示例。上述所有成果可为嵌入式系统设计者提供指导。关于第二点，所选算法已在相对于更高级中央处理器而言硬件资源极为有限的ARM Cortex-M4F微控制器上实现。我们开展了多项实验，以确定在速度、精度和能耗方面性能最优的算法。
关键词 ：奇异值分解 (SVD)；矩阵分解；嵌入式系统；微控制器；ARM Cortex-M；卡尔曼
1. 引言
奇异值分解（SVD）是现代信号处理中基本且最重要的数学工具之一。该方法最初由贝尔特拉米和约旦[1]提出，随后由奥通[2],、埃卡特和杨[3]相继推广。自此以来，SVD已成功应用于大量不同的领域，例如生物医学信号处理[4–10],、图像处理[11–19],、卡尔曼滤波[20–22],、阵列信号处理[23],、动态网络[24],、语音处理[25],、同步定位与地图构建（SLAM）系统[26],以及可变数字滤波器设计[27],等。
在过去的几十年中，人们开发了多种算法来求解奇异值分解[28],，其中戈卢布‐赖施[29],、德梅尔‐卡汉[30],、雅可比旋转[31],、单边雅可比旋转[32],和分治法[33]算法因其在常见应用中能保证良好的性能而被广泛使用。关于奇异值分解及其计算方法的概述可参见[34,35]。
近年来，嵌入式系统中的新兴应用引起了广泛关注，例如多输入多输出（MIMO）系统[36,37],数据分析[38–41],信号稀疏表示[42–47],这些应用要求高效的SVD算法。由于SVD算法归结为求解特征值问题，计算成本较高，因此已提出专用硬件解决方案[48–56]以及并行实现方法[57,58]来克服这一瓶颈。
例如，诸如线性代数软件包 (LAPACK) 和基础线性代数子程序 (BLAS)（在其多种实现中）之类的优化库代表了计算机级CPU上的前沿解决方案，它们采用多种优化方法以高效利用硬件资源。一项重要的优化是重新组织操作流程和数据存储方式，以更好地利用处理器的缓存内存，从而加快从较慢的随机存取存储器（RAM）读取和写入大量数据的速度。
另一项改进来自并行性，即将问题分解为多个较小的子集，并利用现代多核CPU的并行执行能力。单指令多数据 (SIMD) 技术提供了另一种形式的并行性，它通过单条指令对多个相同类型和尺寸的数据元素并行执行相同的操作。例如，一个通常可对两个32位值进行加法运算的CPU，能够在相同时间内并行执行四组8位值的加法运算。因此，SIMD 能更充分地利用微处理器中的可用资源，因为在微处理器上处理大量数据时，部分计算单元处于未使用状态，但仍持续消耗电力。
值得注意的是，所有SVD算法都表现出大量的数据级并行性，因此非常适合通过 SIMD加速和多核并行性来实现。
然而，微控制器中并不存在上述引用的硬件资源（如缓存内存、并行执行单元），因此所有相关优化技术在这一类嵌入式系统上均无用武之地。关于SIMD，实际上 Cortex-M4及更新的嵌入式架构提供了一种初级形式的SIMD指令，但这些指令仅限于对 8位或16位整数值进行操作，因此不适用于需要更高精度浮点数值以保证稳定性的复杂计算。
关于更复杂的架构，文献中已有大量关于如何在现场可编程门阵列 (FPGA) 和其他可编程硬件上实现数学算法的研究[53–56]。这些系统通常较为复杂且昂贵，用于高端应用，并基于大规模并行性的完全不同的计算模型。相反，本文的研究范围集中在经济实惠、低功耗低资源嵌入式系统，这类系统通常基于微控制器，广泛应用于分布式场景，例如物联网概念中。
许多嵌入式系统可用内存非常有限，尤其是数据内存，并且可能不采用底层操作系统以进一步减少资源使用。因此，现有的针对计算机级中央处理器的优化库在嵌入式系统中难以直接使用，也不易移植。由于这一原因，近年来对能够在嵌入式处理器上实现SVD算法的特定加速技术的需求日益增长[59]。
关于SVD实现，存在大量不同的算法，因此从嵌入式系统设计者的角度来看，比较在高性能嵌入式处理器中实现的不同算法所达到的性能具有重要意义。
本文旨在介绍一种在ARM Cortex-M（Arm有限公司，英国剑桥）上的SVD算法的高效实现。更具体地说，由于SVD算法是嵌入式系统中若干新兴应用的核心数学工具，且据我们所知，目前在文献中尚无针对裸机微控制器（即不带操作系统的微控制器）的该算法实现，因此开发适用于这些设备、特别是ARM Cortex-M4F架构[60]的库至关重要。为实现这一目标，我们进行了如下工作：
(i) 对最常见的SVD算法进行全面阐述，以更紧凑和统一的方式编写文献中已有的公式，适合作为综合性参考；
(ii) 在ARM Cortex-M4F（Arm有限公司，英国剑桥）微控制器上实现这些算法，以开发适用于无操作系统的嵌入式系统的库。如前所述，LAPACK和BLAS在计算机级CPU上代表了最先进的解决方案，但由于嵌入式系统可用内存极为有限且缺乏底层操作系统，这些优化库无法直接使用或轻松移植到此类设备上。为此，我们决定实现五种常见的SVD算法（Golub–Reinsch、Demmel–Kahan、Jacobi旋转、单边Jacobi旋转和分治法），以便使其适应资源受限的嵌入式系统并在无操作系统环境下运行。
（iii）寻找最适合裸机嵌入式系统低资源特性的最佳实现方案。为了找出在嵌入式系统中性能最优的算法，本文报告了针对多种实现指标（包括速度、精度和能耗）的性能比较研究。
（iv）通过在惯性测量单元（IMU）的卡尔曼滤波中使用所提出的针对ARM Cortex-M4F优化的库，展示SVD在嵌入式系统中的实际实时应用。该示例还证明了在此类系统上实现SVD的优势，表明相较于传统卡尔曼方法，应用SVD算法可提高卡尔曼滤波的数值精度。
因此，对于有兴趣开发可在低资源微控制器上实现的新算法的研究人员和从业者来说，这篇论文是一个很好的起点。
本文的结构如下。第2节主要聚焦于五种最具代表性的SVD算法，并对它们进行了全面的论述，附录A总结了与奇异值分解变换相关的矩阵代数基本概念。第3节报告了对五种算法性能的比较研究。第4节展示了SVD在IMU数据的卡尔曼滤波中的应用，相较于传统算法实现显示出重要改进。最后，结论部分总结了本项工作。
2. 奇异值分解的算法
设A为一个实(m ×n)矩阵，且m ≥n。已知该分解
$$ A= UΣVT \tag{1} $$
其中
$$ UTU= VTV= VVT= I, Σ= diag(σ1,…, σn) \tag{2} $$
存在[35]。矩阵U由n个对应于AAT的前n个最大特征值的标准正交特征向量组成，而矩阵 V由ATA的标准正交特征向量组成。 Σ的对角元素是AAT特征值的非负平方根，称为奇异值。假设
$$ σ1 ≥ σ2 ≥… ≥ σn ≥ 0 \tag{3} $$
因此，如果 rank(A)=M，则得到 σM+1= σM+2=…= σn= 0。分解（1）称为矩阵 A 的奇异值分解（SVD）。
2.1. 戈卢布-赖恩施算法
该方法由 G. H. Golub 和 C. 赖因施 [29], 提出，直接作用于矩阵 A，从而避免了因计算 ATA 而引起的不必要的数值精度损失。该算法可分为以下连续步骤：
（i）豪斯霍尔德双向对角化（详见附录A.2细节）；（ii）带位移的隐式QR方法（详见附录 A.6细节）；该算法的伪代码见算法 1。
算法1 Golub–Reinsch
要求: A ∈ Rm×n(m ≥ n), ε一个很小的单元舍入误差倍数
使用算法A1计算双对角化。
$$ \left[\begin{array}{c} B \ 0 \end{array}\right] \leftarrow (U_1… U_n)^TA(V_1… V_{n-2}) $$
重复
对于 i= 1，. . . ，n −1执行
•如果 |bi,i+1| ≤ ε( |bii| + |bi+1,i+1|)，则将 bi,i+1 设为零
•找到最大的 q 和最小的 p，使得若 B =
$$ \left[\begin{array}{ccc} B_{11} & 0 & 0 \ 0 & B_{22} & 0 \ 0 & 0 & B_{33} \end{array}\right] $$
Bp Bn−p −q Bq p n−p −q q，则B33为对角矩阵且B22 具有非零超对角线。如果 q = n，则 停止 结束 如果 B22 中有任何对角元素为零，则将同一行中的超对角线元素置零，否则应用 算法A3 结束 结束
2.2. 德梅尔-卡汉算法
该算法的结构基于前面所述的戈卢布‐赖恩施算法。然而，德梅尔‐卡汉算法[30]可以对小奇异值实现更高的精度。该算法包括以下主要连续步骤：
(i) 豪斯霍尔德双向对角化（详见附录A.2获取详细信息）；(ii) 带零位移的QR迭代（详见附录A.7获取详细信息）； 该算法的描述如算法2所示。
算法2 德梅尔‐卡汉
要求： A ∈ Rm×n(m ≥ n) ∈单元舍入误差的一个小倍数
使用算法A1进行双对角化计算。
$$ \left[\begin{array}{c} B \ 0 \end{array}\right] \leftarrow (U_1… U_n)^TA(V_1… V_{n-2}) $$
重复
对于 i= 1: n −1 执行
如果满足相对收敛准则，则将bi,i+1设为零
•找到最大的q和最小的p，使得如果B
$$ \left[\begin{array}{ccc} B_{11} & 0 & 0 \ 0 & B_{22} & 0 \ 0 & 0 & B_{33} \end{array}\right] $$
Bp Bn‐p‐q Bq p n‐p‐q q，则B33为对角矩阵，且B22具有非零的超对角线。若q= n则停止结束
若 B中 22的任意对角线元素为零则将同行的超对角线元素置零否则应用隐式零位移QR算法结束
2.3. 雅可比旋转算法
在这种情况下，给定实对称矩阵 A ∈ R n×n，该算法[31]旨在通过变换获得对角矩阵 B ∈ R n×n
$$ B= J^TAJ \tag{4} $$
其中J表示一系列旋转矩阵。特别是对于第k次旋转操作扫掠可以被重写为
$$ A_{k+1}= J_k^T A_k J_k, A_0= A \tag{5} $$
其中 Jk= J(p,q, θ)是雅可比旋转矩阵，通过角度 θ 旋转矩阵 Ak 的第 p行和第 q行以及第p列和第 q列，使得其 (p,q)和 (q, p)元素为零。在每次迭代步骤中，p和 q的值都被适当选择。参考对应于 p,q列的子矩阵，我们有
$$ \left[\begin{array}{cc} b_{pp} & b_{pq} \ b_{qp} & b_{qq} \end{array}\right]= \left[\begin{array}{cc} c & s \ -s & c \end{array}\right]^T \left[\begin{array}{cc} a_{pp} & a_{pq} \ a_{qp} & a_{qq} \end{array}\right] \left[\begin{array}{cc} c & s \ -s & c \end{array}\right] \tag{6} $$
该算法的关键是确定旋转系数c和s，使得非对角矩阵项bpq和bqp被置为零。
该算法在算法3中报告。
算法3 雅可比旋转
要求：A ∈ Rn×n对称
B ← A ∈ Rn×n
重复
对于 i= 1: n −1 执行
对于 j= i+ 1: n 执行
计算旋转系数 s、c，使得
$$ \left[\begin{array}{cc} b_{ii} & 0 \ 0 & b_{jj} \end{array}\right]= \left[\begin{array}{cc} c & s \ -s & c \end{array}\right]^T \left[\begin{array}{cc} a_{ii} & a_{ij} \ a_{ji} & a_{jj} \end{array}\right] \left[\begin{array}{cc} c & s \ -s & c \end{array}\right] \tag{7} $$
if off(A) < ε，其中 off(A) = √∑i≠j a²ij 和 ε是单位的一个较小倍数舍入 then STOP 结束如果 结束循环 结束循环
2.4. 单边雅可比旋转算法
该算法的主要思想[32]是通过角度θ旋转i列和j列的A，使它们彼此正交。这样， (i, j)元素的ATA被隐式置零，从而得到第i列和第j列的标量积。
设J(i, j, θ)为吉文斯矩阵，当其作用于矩阵A时，可得到
$$ B=(b:1 · · · b:i · · · b:n)= AJ=(a:1 · · · a:i · · · a:n) \left[\begin{array}{cccccc} 1 & \cdots & 0 & \cdots & 0 \ & \ddots & & & \ 0 & \cdots & c & s & \cdots & 0 \ & & -s & c & & \ 0 & \cdots & 0 & \cdots & 1 \end{array}\right] \tag{8} $$
其中i、j列由B给出
$$ b:i= c a:i −s a:j \tag{9} $$
$$ b:j = s a:i+ c a:j \tag{10} $$
该算法的关键在于确定旋转系数c和s，使得乘积B=AJ(i, j, θ)的元素(B^TB)ij和(B^TB)ji被置为零。
该算法的伪代码如算法4所示。
算法4 单边雅可比旋转
要求: A ∈ Rn×n对称
B ← A ∈ Rn×n
重复
对于 i= 1: n −1 执行
对于 j= i+ 1: n 执行
计算旋转系数 s、c，使
$$ \left[\begin{array}{cc} b_{ii} & 0 \ 0 & b_{jj} \end{array}\right]= \left[\begin{array}{cc} a_{ii} & a_{ij} \ a_{ji} & a_{jj} \end{array}\right] \left[\begin{array}{cc} c & s \ -s & c \end{array}\right] \tag{10} $$
如果关闭(A) < ε，其中关闭(A) = √ ∑i≠j a²ij 和 ε 是单位的一个小倍数舍入 then STOP end if end for end for
2.5. 分治算法
对于一个方的对称矩阵 A，奇异值与特征值之间，以及奇异向量与特征向量之间存在关系。事实上，由于 A 可以被对角化，我们可以写出
$$ A= QΛQ^T= A sign(Λ)|Λ|Q^T \tag{11} $$
其中 Λ是特征值对角矩阵，Q是一个酉矩阵（在实数情况下为正交矩阵）。由于我们可以假设 |Λ| 的元素按递减顺序排列，因此对角化（11）对应一个奇异值分解，使得 U = Qsign(Λ)， Σ= |Λ| 且VT= QT。换句话说，奇异值是特征值的绝对值，奇异向量是特征向量（具有相同的范数和方向，但方向不一定相同）。对于非对称矩阵，奇异值和奇异向量与特征值和特征向量没有直接关系，而是与对称矩阵 ATA ∈ R M×M和 AAT ∈ RN×N的特征值和特征向量存在严格的关系。事实上，可以很容易地证明
$$ A^TA= VΣ^2V^T \tag{12a} $$
$$ AA^T= UΣ^2U^T \tag{12b} $$
其中 A被分解为 A= UΣVT。由（12）可得，A的奇异值是矩阵 ATA和 AAT的特征值（后者中等于零的除外）。此外，右奇异向量和左奇异向量分别是 ATA和AAT的特征向量，它们可能相差一个符号（注意只要符号正确，矩阵 A就可以被正确重构）。因此，基于上述分析，奇异值分解可归结为对称矩阵的特征值问题。通常，求解此类问题的方法是迭代方法，包括两个阶段：
(i) 在第一阶段，矩阵 A 被变换为一个矩阵 B，其结构使得特征值和特征向量的计算更加容易。这里假设的一个典型选择是三对角形式。
(ii) 在第二阶段，应用迭代方法来确定特征值和特征向量。
参考第二阶段，分治算法旨在将复杂问题简化为单一问题[33]。该算法适用于维度为 T的 N × N 的三对角且对称的矩阵。
该算法的伪代码在算法5中给出。
算法5 分治法
要求: T ∈ Rn×n三对角的对称
Divide T as
$$ T= \left[\begin{array}{cc} T_1 & 0 \ 0 & T_2 \end{array}\right]+ \rho uu^T $$
对于 i= 1,2执行
计算 Λi,Qi特征值/特征向量，方法如下：
如果 Ti足够小则
直接计算特征值/特征向量
else
递归应用分治算法到 Ti
结束 if
结束循环
使用分解后的 T1, T2 计算 D+ ρvvT，其中
$$ D= \left[\begin{array}{cc} \Lambda_1 & 0 \ 0 & \Lambda_2 \end{array}\right], v= \left[\begin{array}{cc} Q_1^T & 0 \ 0 & Q_2^T \end{array}\right]u \tag{13} $$
通过Li的算法求解 secular 方程 D+ ρvvT= QΛQT 来计算特征值 Λ
Compute eigenvectors as
$$ \left[\begin{array}{cc} Q_1 & 0 \ 0 & Q_2 \end{array}\right]Q $$
3. 实验结果
上一节所述的五种算法已在台式计算机的MATLAB以及Cortex‐M4F微控制器上实现。
这些算法已使用多个矩阵A ∈ R m×n进行了测试。为了确保算法在最普遍的情况下正常工作，选择了若干组尺寸递增的矩阵，每组矩阵具有不同的行列比。考虑到微控制器的内存有限，测试了三组矩阵，其行列比分别为m/n= 1、43和2。进一步增加该比例将导致矩阵无法适应微控制器的有限内存，除非显著降低其秩（内存占用不仅需考虑输入和输出矩阵，还需考虑各种算法所需的中间矩阵和向量）。实际上，优先保证了各组矩阵之间具有可比较的秩。由于RAM内存容量较小，也限制了矩阵的绝对大小，因此与MATLAB中的测试相比，在微控制器上的实验使用了尺寸更小的矩阵。
最大的矩阵A已随机生成，所使用的其他所有矩阵均为完整矩阵A的左上部分。
对于需要对称矩阵的算法，输入矩阵将按如下方式转换为对称形式：首先通过应用算法A1将其转换为双对角形式；然后，设B为双对角矩阵的上三角部分，对BBT应用SVD算法。这样，原矩阵的奇异值即为后者三对角矩阵特征值的平方根。
3.1. MATLAB 实现
这些算法最初在配备 32 GB 随机存取存储器的 6核、12线程英特尔 i7 CPU（英特尔，美国加利福尼亚州圣克拉拉）上使用 MATLAB 实现，以测试其正确性。
表1报告了相同测试中部分样例的精度，其中算法在MATLAB中实现，并以MATLAB内置的SVD函数为基准。所报告的误差是通过 MATLAB内置函数与给定例程计算出的奇异值匹配对之间的相对误差的平均值。
如您所见，除德梅尔‐卡汉算法外，所有算法都能确保与MATLAB内置的奇异值分解函数相同的精度。这是可以预期的，因为德梅尔‐卡汉算法专门设计用于处理非常小的奇异值，而我们在实验中使用的随机生成矩阵的奇异值接近于1。
表1. MATLAB中实现的SVD算法与MATLAB内置函数的平均误差。
| Size | 戈卢布‐赖因施 | 德梅尔‐卡汉 | 雅可比旋转 | 单边雅可比 | 分治法 |
|---|---|---|---|---|---|
| 32 × 24 | 0 | 2.0 × 10⁻⁸ | 0 | 0 | 0 |
| 48 × 36 | 0 | 2.7 × 10⁻⁸ | 0 | 0 | 0 |
| 96 × 72 | 0 | 5.4 × 10⁻⁸ | 0 | 0 | 0 |
| 128 × 96 | 0 | 7.9 × 10⁻⁸ | 0 | 0 | 0 |
| 160 × 120 | 0 | 8.6 × 10⁻⁸ | 0 | 0 | 0 |
| 200 × 150 | 0 | 1.9 × 10⁻⁷ | 0 | 0 | 0 |
3.2. Cortex-M4F 实现
在位于瑞士日内瓦的意法半导体（STMicroelectronics）生产的STM32F429ZI微控制器（https://www.st.com/en/microcontrollers-microprocessors/stm32f429zi.html）上测试了不同算法的Cortex‐M4F实现，该微控制器安装在32F429IDISCOVERY评估板上（ https://www.st.com/en/evaluation-tools/32f429idiscovery.html）。
STM32F429ZI 微控制器基于 ARM 32位 Cortex-M4F [60] CPU，主频为 180兆赫，具有 2兆字节的闪存用于代码和只读数据，以及 256千字节的随机存取存储器。此外，它还具有多个与此工作无关的硬件外设。
Cortex-M4F 内核的一个主要特点是内置了 32位 硬件浮点运算单元（FPU），其名称中额外的“F”即表示此含义。对于在浮点域中进行的各种高计算量工作而言，FPU 至关重要，本文中进行的 SVD 实验正是这种情况。Cortex-M4F 的 FPU 仅支持 32位（ https://developer.arm.com/docs/ddi0439/latest/floating‐point‐unit/about‐the‐fpu），因此算法采用单精度（32位）值实现。若在没有 FPU 或需要比硬件支持更高精度的情况下实现此类算法，则必须使用软件浮点数学库，这在资源本已受限的系统中将是不可行的。
算法使用C语言实现。未使用特定的开发环境，代码在GNU‐Linux机器上使用GCC软件套件针对ARM进行编译，采用自定义的makefile，并借助意法半导体（瑞士日内瓦）提供的STM32F4xx标准外设驱动，这是一套由ST为其特定微控制器提供的库，涵盖了从底层初始化到硬件外设使用的所有硬件管理方面。固件为“裸机”类型，因此未添加任何实时操作系统（RTOS）或其他中间件。
硬件系统没有提出任何特殊设计要求，仅使用了32F429IDISCOVERY开发板（意法半导体，瑞士日内瓦）已提供的功能。该设备在其最高速度180兆赫下运行。该开发板还集成了对微控制器进行编程与调试所需的所有硬件，即ST‐LINK/V2接口（意法半导体，瑞士日内瓦），可通过USB与计算机连接。在计算机端，通过使用OpenOCD（ http://openocd.org）实现与该接口的通信，OpenOCD是一款用于ARM及其他系统的免费调试与编程软件。最终，OpenOCD作为GCC套件中标准调试器GDB的服务器，在需要时用于将代码下载到设备并检查其内存以获取测试结果。
关于输入和输出，程序的只读数据（例如从中生成较小矩阵的较大矩阵，或用于计算精度的奇异值参考向量）存储在比RAM更充足的程序内存（闪存）中。一旦程序运行一系列测试，可以通过调试器在合适的时间点中断程序来检查数值输出。单个例程的时序由软件本身使用Cortex-M内核内置的SysTick定时器进行计算。
除了编译器执行的优化之外，还特别注重优化关键的数学例程，同时尽量减少程序内存和所需RAM的使用，以在系统受限的资源范围内尽可能加快计算速度。
一种初步的优化方法是选择内存中数据的最佳排列方式。二维矩阵可以以两种不同的方式存储在随机存取存储器中（即字节的线性数组）：行优先或列优先顺序，分别表示按行或按列依次存储数据（参见图2）。前者是在C语言中或任何遵循数学惯例（将行作为第一索引）的情况下最常用的方式。这种方式在带有缓存内存的中央处理器中具有更为关键的影响，因为缓存会从随机存取存储器加载顺序数据，从而使得访问顺序数据的速度远高于稀疏数据。如前所述，微控制器没有缓存内存，因此这种情况并不直接适用；然而，由于具有自动递增功能的加载/存储汇编指令的存在，在顺序访问数据时仍具有不可忽视的优势，这类指令能够在单个执行周期内读取或写入内存数据并自动递增地址寄存器。
作为行优先或列优先数据存储影响的一个定量示例，我们可以考虑单边雅可比旋转算法，该算法与其他算法的不同之处在于它仅按列访问输入矩阵。表2显示了针对最大矩阵集合的算法在绝对时间和速度提升百分比方面的不同时间。如前所述，速度的提升虽然可观，但幅度很小。最后一列表示在矩阵以相反方式存储的情况下转置矩阵所需的时间，该时间大约比获得的时间改进小一个数量级。此外，SVD 的输入矩阵通常来自先前的计算，因此可以一开始就以更合适的方式生成。
表2. 单边雅可比旋转在内存中数据排列方式不同时的时序差异。
| 矩阵大小 | 行优先 | 列主序 | 速度提升 | 转置 |
|---|---|---|---|---|
| 48 × 24 | 19.6 毫秒 | 19.1 毫秒 | 2.6% | 0.07 毫秒 |
| 72 × 36 | 73.1 毫秒 | 71.9 毫秒 | 1.7% | 0.15 毫秒 |
| 96 × 48 | 192 毫秒 | 189.7 毫秒 | 1.2% | 0.25 毫秒 |
| 120 × 60 | 370.6 毫秒 | 366.8 毫秒 | 1% | 0.39 毫秒 |
| 144 × 72 | 634 毫秒 | 628.6 毫秒 | 0.9% | 0.55 毫秒 |
除了以方便的顺序访问数据（这会带来适度的速度提升）以及缺乏可进一步利用的硬件资源外，其他优化必须通过尽可能减少不必要的计算来实现。一个典型的例子是矩阵乘法，这是计算成本最高的操作之一矩阵代数中的运算。一般来说，如果C= A · B，则C的通用元素由下式给出
$$ c_{i,j}= \sum_{k=1}^{n} a_{i,k}b_{k,j} \tag{14} $$
其中n是A的列数。以这种方式计算C的所有元素需要一个三重嵌套循环，计算成本非常高，特别是对于大型矩阵而言。
通过观察特定矩阵的性质，可以减少矩阵乘法所需的运算次数。例如，公开算法中的一个重复操作要求对方阵与其转置进行矩阵乘法，如C= A · A′。在这种情况下，(14) 变为
$$ c_{i,j}= \sum_{k=1}^{n} a_{i,k}a_{j,k}. \tag{15} $$
首先我们可以注意到，在内层循环中A总是按行遍历，因此如果矩阵以行主序存储，则可以始终以最方便的顺序读取数据。更重要的是，可以明显看出A · A′是一个对称矩阵，因此我们实际上只需计算其大约一半的元素，从而显著减少运算次数（图 3）。
在C= A · A′的情况下，可以进一步减少运算次数，其中A是上对角矩阵，这是给定算法中的另一个常见情况。由于ai,j ≠ 0仅当j= i或j= i+ 1时成立，由(15)可得，ci,j ≠ 0仅当j= i −1，i或i+ 1时成立（实际上结果矩阵是三对角矩阵），并且求和中唯一的非零项是那些ai,k和aj,k都为非零的项。最终公式为：
$$ c_{i,j} = \begin{cases} a_{i,i}a_{j,i}+ a_{i,i+1}a_{j,i+1}, & j= i \ a_{i,i+1}a_{j,i+1}, & j= i+ 1 \ a_{i,i}a_{j,i}, & j= i−1 \ 0, & \text{otherwise} \end{cases} \tag{16} $$
其中 (i+ 1) th元素可能不存在。由于该矩阵也是对称的，因此前述情况的约简同样适用 （图 4）。
另一个特殊的例子是矩阵与吉文斯矩阵的乘法，形式如（A14），用于执行吉文斯旋转。为避免与下标混淆，我们称 G为吉文斯矩阵，并为简化起见仅限于左乘的情况，如 C= G· A。如果初始时设 C= A，则根据 G的定义可知，只有 p和 q行的 C会受到乘法的影响。此外，在计算中仅使用矩阵 A某一列中的第 p和 q个元素（图5）。因此，只需更新 C中位于 p和 q行的元素，其值为：
$$ c_{p,j}= g_{p,p}a_{p,j}+ g_{p,q}a_{q,j}, c_{q,j}= g_{q,p}a_{p,j}+ g_{q,q}a_{q,j} \tag{17} $$
对于每一列 j。右乘也有类似的公式。
先前计算的复杂度（对应于输入矩阵不同特殊情况下的矩阵乘积）可以通过输入矩阵尺寸相关的标量乘法次数进行比较。结果如表3所示，并附有特定软件例程的代码尺寸。
表3。 不同类型的 n ×n输入矩阵在矩阵乘积中的标量乘法次数。
| 矩阵乘积 | 标量乘积 | 代码大小（字节） |
|---|---|---|
| C= A · B，通用 | n³ | 148 |
| C = A · A′ | n³/2+ n²/2 | 210 |
| C= A · A′，A 上对角 | 3n −1 | 304 |
| C = G· A,吉文斯旋转 | 4n | 108 |
此类优化对计算速度的影响可以进行测量，特别是来自（15）和（16）的优化以及对称性特性带来的优化（吉文斯旋转实际上从未以矩阵乘法实现，因为其矩阵明确构造为仅修改结果中的少数元素）。表4展示了在使用两种算法时，针对最大矩阵集获得的速度提升，这些优化在这些算法中最为重要。
表4. 优化矩阵乘法算法带来的速度提升。
| 矩阵大小 |
表5. Cortex-M4F上SVD算法的时间。
| Size | 戈卢布- 赖因施 | 德卡梅汉尔-雅可比旋转 | 单边J. | 征分服治 |
|---|---|---|---|---|
| 24 × 24 | 17 ms | 22 毫秒 | 28 毫秒 | 13 毫秒 |
| 36 × 36 | 78 毫秒 | 133 毫秒 | 112 毫秒 | 45 毫秒 |
| 48 × 48 | 229 毫秒 | 441 毫秒 | 310 毫秒 | 114 毫秒 |
| 60 × 60 | 681 毫秒 | 1100 毫秒 | 693 毫秒 | 218 毫秒 |
| 72 × 72 | 1092 毫秒 | 1232 毫秒 | 1396 毫秒 | 408 毫秒 |
| 32 × 24 | 28 毫秒 | 34 毫秒 | 39 毫秒 | 16 毫秒 |
| 48 × 36 | 130 毫秒 | 186 毫秒 | 164 毫秒 | 57 毫秒 |
| 64 × 48 | 388 毫秒 | 463 毫秒 | 482 毫秒 | 131 毫秒 |
| 80 × 60 | 939 毫秒 | 1272 毫秒 | 1100 毫秒 | 252 毫秒 |
| 96 × 72 | 2963 毫秒 | 2878 毫秒 | 2182 毫秒 | 430 毫秒 |
| 48 × 24 | 58 毫秒 | 64 毫秒 | 69 毫秒 | 19 毫秒 |
| 72 × 36 | 283 毫秒 | 299 毫秒 | 316 毫秒 | 72 毫秒 |
| 96 × 48 | 1115 毫秒 | 1434 毫秒 | 970 毫秒 | 190 毫秒 |
| 120 × 60 | 2429 毫秒 | 2280 毫秒 | 2299 毫秒 | 367 毫秒 |
| 144 × 72 | 4427 毫秒 | 4669 毫秒 | 4672 毫秒 | 628 毫秒 |
表6. SVD算法在Cortex-M4F与MATLAB内置函数上的平均误差。
| Size | 戈卢布- 赖因施 | 德梅尔- 卡汉雅可比旋转 | 单边 J. | 分治 征服 |
|---|---|---|---|---|
| 24 × 24 | 3.8 × 10⁻⁷ | 7.8 × 10⁻⁷ | 1.0 × 10⁻⁶ | 1.9 × 10⁻⁷ |
| 36 × 36 | 4.8 × 10⁻⁷ | 2.3 × 10⁻⁶ | 1.7 × 10⁻⁶ | 3.5 × 10⁻⁷ |
| 48 × 48 | 4.7 × 10⁻⁷ | 7.1 × 10⁻⁶ | 3.0 × 10⁻⁶ | 2.4 × 10⁻⁷ |
| 60 × 60 | 5.3 × 10⁻⁷ | 4.8 × 10⁻⁶ | 7.3 × 10⁻⁷ | 3.0 × 10⁻⁷ |
| 72 × 72 | 4.6 × 10⁻⁷ | 2.6 × 10⁻⁶ | 1.7 × 10⁻⁶ | 3.4 × 10⁻⁷ |
| 32 × 24 | 2.7 × 10⁻⁷ | 5.7 × 10⁻⁷ | 2.2 × 10⁻⁷ | 1.7 × 10⁻⁷ |
| 48 × 36 | 3.6 × 10⁻⁷ | 2.4 × 10⁻⁶ | 4.0 × 10⁻⁷ | 1.7 × 10⁻⁷ |
| 64 × 48 | 2.9 × 10⁻⁷ | 2.7 × 10⁻⁶ | 2.9 × 10⁻⁷ | 1.7 × 10⁻⁷ |
| 80 × 60 | 4.1 × 10⁻⁷ | 5.9 × 10⁻⁶ | 5.4 × 10⁻⁷ | 2.4 × 10⁻⁷ |
| 96 × 72 | 1.6 × 10⁻⁶ | 6.5 × 10⁻⁶ | 6.0 × 10⁻⁷ | 2.7 × 10⁻⁷ |
| 48 × 24 | 3.0 × 10⁻⁷ | 1.7 × 10⁻⁶ | 1.6 × 10⁻⁷ | 1.7 × 10⁻⁷ |
| 72 × 36 | 2.7 × 10⁻⁷ | 6.5 × 10⁻⁷ | 3.7 × 10⁻⁷ | 1.5 × 10⁻⁷ |
| 96 × 48 | 1.0 × 10⁻⁶ | 8.7 × 10⁻⁶ | 4.1 × 10⁻⁷ | 1.8 × 10⁻⁷ |
| 120 × 60 | 5.9 × 10⁻⁷ | 3.9 × 10⁻⁶ | 4.8 × 10⁻⁷ | 2.0 × 10⁻⁷ |
| 144 × 72 | 4.8 × 10⁻⁷ | 5.3 × 10⁻⁶ | 5.8 × 10⁻⁷ | 3.1 × 10⁻⁷ |
表7. Cortex-M4F上SVD算法的能耗。
| Size | 戈卢布- 赖因施 | 德梅尔- 卡汉 | 雅可比旋转 | 单边 J. | 分治 征服 |
|---|---|---|---|---|---|
| 32 × 24 | 5.7 毫焦 | 7.1 毫焦 | 10.9 毫焦 | 1.2 毫焦 | 5.8 毫焦 |
| 48 × 36 | 26.7 毫焦 | 30.6 毫焦 | 38.3 毫焦 | 9.0 毫焦 | 26.1 毫焦 |
| 64 × 48 | 79.7 毫焦 | 87.8 毫焦 | 113.6 毫焦 | 26.5 毫焦 | 78.8 毫焦 |
| 80 × 60 | 193.5 毫焦 | 217.5 毫焦 | 250.7 毫焦 | 50.9 毫焦 | 187.2 毫焦 |
| 96 × 72 | 594.6 毫焦 | 1150.9 毫焦 | 487.8 毫焦 | 86.7 毫焦 | 382.1 毫焦 |
4. 应用：惯性系统数据的卡尔曼滤波
作为嵌入式系统的一个实际应用，并为了证明在此类系统上实现SVD的优势，我们将展示如何利用SVD提高应用于惯性测量单元（IMU）的卡尔曼滤波的数值精度。IMU中包含的传感器通常为加速度计和陀螺仪。通过对这些传感器的数据进行适当的滤波和数据融合技术处理，可以计算出运动物体的空间姿态（即其姿态），也就是物体相对于空间轴的相对角度。
加速度计存在信号上叠加高频噪声的问题，且当对物体施加力时还会引入额外的加速度项。陀螺仪测量的是围绕各轴的角速度变化率，必须通过积分才能得到角位置；由于陀螺仪的测量受到非恒定且不可预测的偏差影响，导致所获得的数据会出现明显的漂移，且漂移随时间逐渐增加。
卡尔曼滤波器[61]是一种常用的技术，用于根据随时间变化且受统计噪声影响的离散测量值来估计此类系统（即沿选定轴的当前角度）的状态。为此，必须将系统动态表示为状态空间方程的形式：
$$ x_t= Ax_{t-1}+ Bu_{t-1}+ w_{t-1} \tag{18} $$
$$ y_t= Cx_t+ v_t \tag{19} $$
其中，$x_t$ 是待估计的未知系统状态向量，向量 $u_t$ 是系统的控制输入，$y_t$ 是测量向量。$w$ 和 $v$ 分别是过程和测量数据上的加性噪声，均假设为零均值高斯过程：$w_t ∼ N(0, Q_t)$，$v_t ∼ N(0,R_t)$。初始状态条件 $x_0$ 也是一个随机变量：$x_0 ∼ N(0,\Pi_0)$。
[22]表明，传统的卡尔曼滤波实现（如基于乔列斯基分解的平方根算法）在处理良态问题时是有效的，但由于这些算法在多处依赖矩阵求逆，当方程为病态时，可能会因舍入误差而失效。另一方面，[22]展示了基于奇异值分解的卡尔曼滤波变体在病态情况下表现更好。
在IMU示例中，病态问题的一种情况是存在多组相同的测量值。为了更好地展示一个实际案例，我们报告了通过以下方式获得的结果
一块 X-NUCLEO-IKS01A2（https://www.st.com/en/ecosystems/x-nucleo-iks01a2.html） 板，这是意法半导体推出的一款传感器板，配备有陀螺仪-加速度计组合传感器和磁力计-加速度计等。该板可与微控制器单元配合使用，在本例中我们使用了 NUCLEO-F411RE（ https://www.st.com/en/evaluation-tools/nucleo-f411re.html），其搭载基于 Cortex-M4F 的微控制器，与之前实验中使用的微控制器类似。
对于该硬件系统，固件使用 STM32CubeIDE 开发平台（版本 1.4.0）进行开发。通过集成的软件包管理器，在 IDE 中安装了额外的软件包 X-CUBE-MEMS1。该软件包支持与 X-NUCLEO-IKS01A2 惯性系统中的特定传感器进行通信。在此情况下，微控制器（ STM32F411RE）的时钟频率为 84 兆赫兹。与之前的实验一样，该固件为“裸机”模式，未添加 RTOS 或其他中间件。本例中，与开发板之间的编程与调试通信完全由 IDE 管理。
如果有两个加速度计可用，利用这两个加速度计的数据来获得更优的状态估计是较为方便的，因为器件本身具有噪声特性。从数学上讲，公式(18)和(19)中的矩阵可以表示为：
$$ A= \left[\begin{array}{cccc} 1 & -\Delta t & 0 & 0 \ 0 & 1 & 0 & 0 \ 0 & 0 & 1 & -\Delta t \ 0 & 0 & 0 & 1 \end{array}\right] B= \left[\begin{array}{cc} \Delta t & 0 \ 0 & 0 \ 0 & \Delta t \ 0 & 0 \end{array}\right] C= \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \ 0 & 0 & 1 & 0 \ 1+ \delta & 0 & 0 & 0 \ 0 & 0 & 1+ \delta & 0 \end{array}\right] \tag{20} $$
以及变量如下：
$$ x= \left[\begin{array}{c} \hat{\varphi} \ \hat{b}
\theta \ \hat{\varphi} \ \hat{b}
\theta \end{array}\right] u=\left[ \dot{\varphi}
G \dot{\theta}_G \right] y= \left[\begin{array}{c} \hat{\varphi}
{A1} \ \hat{\theta}
{A1} \ \hat{\varphi}
{A2} \ \hat{\theta}
{A2} \end{array}\right] \tag{21} $$
此处系统状态包括沿x轴和y轴的估计横滚角和俯仰角（$\hat{\varphi}, \hat{\theta}$），以及这两个角度上的估计陀螺仪偏置（$\hat{b}
\varphi, \hat{b}
\theta$）；系统输入包括来自陀螺仪的角速度（ $\dot{\varphi}_G, \dot{\theta}_G$）； $\Delta t$为采集数据周期；系统测量值由来自两个加速度计数据的估计角度给出（$\hat{\varphi}
{A1}, \hat{\theta}
{A1}, \hat{\varphi}
{A2}, \hat{\theta}_{A2}$）。
$\delta$因子反映了两组加速度计数据在数值上非常相似，从而导致病态问题：由于系统测量变量的冗余，矩阵接近奇异。如[22],所示，当 $\delta$小到与机器精度相当时，可以预期使用传统算法求解卡尔曼滤波方程将受到舍入误差的影响。我们还可以预期，在这种情况下，基于奇异值分解的算法将表现得更好。
板载固件已用于将一批测量数据（陀螺仪、加速度计1、加速度计2）导出到计算机，以便能够在MATLAB中以可重复的方式对数据进行检查。数据的持续时间约为90秒， $\Delta t$ 为10毫秒，通过将设备在三个轴上的若干参考位置（即0度、90度和 −90度）移动而获得。
然后，在MATLAB中使用两种算法（常规算法和基于奇异值分解的算法）计算了卡尔曼滤波器，并针对 $\delta$的不同递减值重复实验，直至接近机器精度的量级。
最后，通过计算所有不同 $\delta$值下估计状态相对于已知方向参考区间的平均绝对误差 （单位：度），来评估其精度。
表 8显示，这两种算法在降至某个 $\delta$值时结果非常接近，但在更低的值下，传统算法失效 （“NaN”表示非数字，是一种指示无效数值的错误状态）。相反，基于奇异值分解的算法在 $\delta$的另外两个数量级范围内，仍能持续工作且差异可忽略，仅在更低的值时才开始退化。
表8。 MATLAB：估计角度的平均绝对误差（φ, θ）（度）。
| 算法 | δ= 10⁻⁸ | δ= 10⁻⁹ | δ= 10⁻¹⁰ | δ= 10⁻¹¹ | δ= 10⁻¹² |
|---|---|---|---|---|---|
| 传统的 | 1.7337 1.1066 | 1.7337 1.1066 | NaN NaN | NaN NaN | NaN NaN |
| 基于奇异值分解 的 | 1.7337 1.1066 | 1.7337 1.1066 | 1.7357 1.1065 | 1.7113 1.3738 | 23.6834 27.0080 |
为了在嵌入式系统有限的能力下进行类似的实验，已在相同的离线数据集上用C程序实现了相同的卡尔曼滤波器算法，使用了前几节讨论的自定义SVD算法，并采用单精度 （32位）浮点数，即与Cortex-M4F架构的浮点硬件相同的精度。为了仅关注传统卡尔曼滤波器与基于奇异值分解的卡尔曼滤波器之间的差异，我们只使用了一种SVD算法，即单边雅可比算法，该算法已被证明在多个标准下是最佳选择（第3节）。其他算法显示出类似的结果。
初步测试表明，在数值精度受限的情况下（如此处所示），传统的卡尔曼算法对初始条件的变化更为敏感。实际上，即使是在单个加速度计的情况下，通过改变所涉及的协方差矩阵的相对值，两种卡尔曼实现之间的差异也是明显的。
作为一个说明性示例，展示了一系列测试，其中 $\Pi_0= I_4$和$Q= I_4$,使用单个加速度计的固定值数据，而 $R$取不同值。结果如表 9所示，显示出与MATLAB情况类似的趋势： 两种算法在某些情况下非常接近，但基于奇异值分解的算法在更广泛的参数范围内表现一致。
表9. 单精度C代码：估计角度的平均绝对误差（度）(φ, θ)。
| 算法 | R= I₂ | R = 0.8I₂ | R = 0.6I₂ | R = 0.4I₂ | R = 0.2I₂ |
|---|---|---|---|---|---|
| 传统的 | 2.7197 1.3642 | 2.7191 1.3820 | 2.7124 1.4480 | NaN NaN | NaN NaN |
| 基于奇异值分解 的 | 2.7186 1.3477 | 2.7173 1.3494 | 2.7159 1.3516 | 2.7143 1.3543 | 2.7124 1.3580 |
关于在Cortex-M4F微控制器上的实际性能，基于奇异值分解的卡尔曼算法所涉及的矩阵比表5中的矩阵要小得多，在双加速度计情况下最大尺寸为8x4。对表5中结果的外推表明计算耗时仅为毫秒的几分之一，因此在大多数应用中，使用基于奇异值分解的变体实现实时计算卡尔曼滤波器是可行的。
5. 结论
本文对五种最常见的SVD算法进行了比较研究，即Golub–Reinsch、Demmel– Kahan、雅可比旋转、单边雅可比旋转和分治法。该比较研究的目的是找出最适合嵌入式系统的算法。所选算法从理论角度进行了研究，并给出了详细的数学描述，为实际实现提供了关键提示。这些算法最初用MATLAB编写以验证其正确性，随后在Cortex-M4F微控制器上实现，以获得针对低资源嵌入式系统的优化结果。对Cortex-M4F实现的比较表明，单边雅可比旋转在速度、精度和能耗方面优于其他算法。此外，分治算法表现出相似的精度结果。最后，将SVD应用于一个示例案例，表明在数据为病态的情况下，它可以提高算法的精度，而传统实现可能会失败。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
