通过示例了解 C++11 中的奇异值分解 (SVD)
5.00/5 (14投票s)
通过示例解释 C++11 中的 SVD
引言
奇异值分解 (Singular Value Decomposition, SVD) 是将实矩阵分解为规范形式的方法。奇异分解是在处理矩阵时一种方便的方法。它展示了矩阵的几何结构，并允许您可视化可用数据。奇异分解用于解决各种问题——从最小二乘法的近似和方程组的求解到图像压缩。同时，利用了奇异分解的不同性质，例如，显示矩阵秩的能力，近似给定秩的矩阵。SVD 允许计算大型矩阵的逆矩阵和伪逆矩阵，使其成为解决回归分析问题的有用工具。
奇异值分解最初由数学家发明并提出，旨在确定给定矩阵经过两个空间各种独立正交变换后是否具有另一种自然的双线性形式。
在1873年和1874年，Eugenio Beltrami 和 Camille Jordan 揭示，将双线性形式表示为矩阵的奇异值，在进行正交替换时，会产生一套完整的双线性形式不变量。James Joseph Sylvester 也在1889年发明并提出了实方矩阵的奇异值分解。
第四位独立建立奇异值分解的数学家是 Autonne，他在1915年通过极分解计算了 SVD。Carl Eckart 和 Gale Young 在1936年首次证明了矩形和复数矩阵的奇异值分解。他们将其视为厄米特矩阵主轴变换的推广。
1907年，Erhard Schmidt 为积分算子（在某些弱技术假设下是紧的）定义了奇异值的模拟；他似乎不熟悉有限矩阵奇异值的并行工作。Émile Picard 在1910年进一步提出了这一理论，他是第一位将计算出的数值称为“奇异值”的数学家。
计算SVD的实用方法可以追溯到Kogbetliantz在1954年、1955年以及Hestenes在1958年，它们与Jacobi特征值算法非常相似，都使用平面旋转或Givens旋转。然而，这些方法后来被Gene Golub和William Kahan在1965年发表的方法所取代，后者使用Householder变换或反射。在1970年，Jacque Golub和Christian Reinsch发表了Golub/Kahan算法的一个变体，该变体至今仍是最常用的。
在本文中，我们将讨论用于计算给定实数或酉矩阵完整 SVD 的简单迭代法和 Jordan-高斯变换法。此外，我们将演示用 C++11 实现 SVD 计算算法的代码，并进行深入讨论。具体来说，我们将提供一个分步教程，以一个积分矩阵的奇异值分解为例，计算完整的 SVD。
背景
计算奇异值分解（SVD）
矩阵A的奇异值分解（SVD）是一种算法，它允许我们将给定的实矩阵或复矩阵A分解为一组奇异值，以及其左奇异向量和右奇异向量。代数上，奇异值分解可以表述为
$\begin{aligned}A=U\ast S\ast V^T\end{aligned}$
其中，\(A\) – 是一个给定的实矩阵或酉矩阵，\(U\) – 左奇异向量的正交矩阵，\(S\) – 是奇异值的对称对角矩阵，\(V^T\) – 分别是右奇异向量的转置正交矩阵。
从上面的公式 (1) 可以看出，给定矩阵的分解 (A) 是左奇异向量的某个正交矩阵 \(U\)、对称对角奇异值矩阵 \(S\) 和右奇异向量的转置正交矩阵 \(V^T\) 的乘积
根据上面列出的公式，\(\forall s_i\)\((s_1\geq\ s_2\geq\ s_3\geq\ldots\geq\ s_{\min(m,n)})\) 的每个奇异值精确地对应一对奇异向量，可以是 \(U_i=\left\{\begin{matrix}u_{1,1}&u_{1.2}&u_{1,3}\\\end{matrix}\ldots\begin{matrix}u_{1,(n-1)}&u_{1,n}\\\end{matrix}\right\}\) 或 \(V_i=\left\{\begin{matrix}v_{1,1}&v_{2,1}&v_{3,1}\\\end{matrix}\ldots\begin{matrix}v_{(m-1),1}&v_{m,1}\\\end{matrix}\right\}\)。形式上，给定矩阵 \(A\) 的 SVD 分解必须满足以下条件，使得
其中 \(A\) – 是一个给定矩阵，\(U_i\) 和 \(V_i\) – 左、右正交奇异向量，\(s_i\) – 奇异值。通过将第二个公式应用于上面列出的第一个公式，我们可以证明 \(i-th\) 向量 \(U_i\) 和 \(V_i\) 实际上是特征向量，并且 \(i-th\) 值 \(s_i \) – 是恰好对应这两个向量的特征值。
步骤 1：计算对称分解矩阵 \(A^TA\)
要计算矩阵 \(A\) 的完整 SVD，我们必须首先计算某个对称实矩阵或酉矩阵的特征值和特征向量，该矩阵是矩阵 \(A^T\) 和 \(A\) 的乘积。这两个矩阵的乘积实际上给我们一个对称矩阵，其特征值 \(\sigma\) 很容易计算。
例如，假设我们给定一个矩阵 \(A\)
让我们找到一个对称分解矩阵 \(A^TA\)
成功计算 \(A^TA\) 后，我们来求这个对称分解矩阵的特征值。显然，分解矩阵将具有以下特征值：\(\sigma_1=15.4310,\ \sigma_2=5.5573,\ \sigma_3=0.0116\)。
步骤 2：计算矩阵 \(A^TA\) 的奇异值
由于分解矩阵 \(A^TA\) 的一组特征值已经计算出来，现在我们来求奇异值对角矩阵 \(S\)。为此，我们必须计算每个特征值的平方根，并将它们放置在矩阵 \(S\) 的对角线上
之后，我们还必须计算逆矩阵 \(S^{-1}\)。正如我们已经发现的，所有奇异值必须以降序排列在矩阵 \(S\) 的对角线上，从第一个最大奇异值 \(s_{max}\) 开始。
为了找到逆对角矩阵 \(S^{-1}\)，我们所要做的就是将 1.0 除以矩阵 \(S\) 对角线上的每个奇异值
例如，假设我们已经计算了步骤1中的那些特征值。现在我们计算奇异值矩阵 \(S\) 及其逆矩阵 \(S^{-1}\)
步骤 3：计算矩阵 \(A^TA\) 的右奇异向量
既然我们已经计算了奇异值的对角矩阵并找到了它的逆，现在我们来获取 \(A^TA\) 分解矩阵的右特征向量，并通过将每个右特征向量的分量除以每个向量的绝对标量长度 \(L\) 来找到矩阵 \(A\) 的右奇异向量 \(V\)
此外，为了计算矩阵 \(A\) 的右奇异向量，我们还必须在执行这些计算之后，找到所得矩阵 \(V\) 的转置。
例如，假设我们已经使用公式 \((A^TA-\sigma I)=0\) 计算了分解矩阵 \(A^TA\) 的右特征向量
这三个右特征向量中的每一个都可以通过找到给定线性方程组的非平凡解来获得，通过使用 Jordan-Gaussian 方法求解，该方法允许我们将给定矩阵转换为简化行阶梯形，如下文所述。
显然，分解矩阵 \(A^TA\) 将具有以下右特征向量 \(V\)
现在我们来找到这些向量的绝对标量长度
最后，我们将这三个向量沿行放置在矩阵 \(V\) 中，并转置所得矩阵
显然，转置矩阵 \(V^T\) 是矩阵 \(A\) 的右奇异向量的正交矩阵。
步骤 4：计算矩阵 \(A^TA\) 的左奇异向量
最后，由于我们已经获得了右奇异值的正交矩阵，我们可以很容易地将给定矩阵 \(A\)、右奇异值的转置矩阵 \(V^T\) 和先前获得的逆对角矩阵 \(S^{-1}\) 的乘积计算出左奇异值 \(U\) 的特定矩阵。这通常通过以下公式完成
应用上面列出的特定变换后，矩阵 \(U\) 的每一行都将包含计算出的左奇异向量。例如，让我们使用上面列出的公式找到矩阵 \(A\) 的左奇异向量正交矩阵 \(U\)
执行所有这些计算的结果，我们将得到矩阵 \(A\) 的以下完整 SVD 分解
在本文的下一段中，我们将具体讨论如何计算给定对称分解矩阵 \(A^TA\) 的特征值和特征向量，并提供执行此类计算的详细示例。
特征值和特征向量（EVD）
正如我们已经讨论过的，给定矩阵 A 的完整 SVD 分解过程主要基于特征值和特征向量的计算。实际上，为了计算给定矩阵 A 的这些特征向量，我们必须首先计算该矩阵的所有特征值。
在线性代数中，有许多方法可以让我们找到小型矩阵的特征值，例如 2x2、3x3 等。例如，为了找到给定矩阵的所有特征值，我们可以计算其行列式。找到矩阵 A 行列式的结果是一个 N 次方方程，其解实际上就是我们需要找到的那些特征值。然而，大多数现有方法用于寻找矩阵 A 特征值并非计算性的，或者不能被公式化为计算机算法。
在这一段中，我们将讨论一种称为简单迭代法的方法，它允许我们找到分解矩阵 \(A^TA\) 的所有特征值，然后找到与每个获得的特征值精确对应的该矩阵的每个特征向量。与上面回顾的其他现有方法不同，以下方法可以很容易方便地被公式化为计算算法。
一般来说，使用简单迭代法的思想如下图所示
步骤 1：计算矩阵 \(A^\ast\) 的最大特征值
我们首先需要做的是找到给定矩阵 \(A^\ast=A^TA\) 的第一个最大特征值。为此，我们必须使用简单迭代法，该方法可以表述如下：
-
\(A_{(mxm)}^\ast=A^TA\) – 是一个分解矩阵，\(R=\left\{1,1,1\ldots,1\right\}\) – 是一个初始单位向量。此外，设 \(i=\bar{1,t}\) – 是一个循环计数器变量，\(\varepsilon=10e-6\) – 是一个精度误差的常数值，\(V^\ast\) – 是一个存储每次迭代结果的向量，\(M_{(mxt)}\) – 是一个结果矩阵，其列是向量 \(V^\ast\) 的值，\(\ \ \sigma\) 和 \(\sigma_{old}\) – 分别是当前迭代和前一次迭代的特征值 (\(\sigma=0,\ \sigma_{old}=\sigma)\))；
-
找到矩阵 \(A^\ast\) 和特定向量 \(R\) 的乘积，并得到向量 \(V^\ast=A^\ast\ast R\)；
-
通过将向量 \(V^\ast\) 的值沿其 \(i-th\) 列放置，将其附加到结果矩阵 \(M\)；
-
通过将结果向量 \(V^\ast\) 的值赋给向量 \(R\)\((R=\ V^\ast)\) 来更新向量 \(R\)；
-
如果这不是第一次迭代 \((i>0)\)，则通过对结果矩阵 \(M\) 的第一行中的 \(i-th\) 和 \((i\ – 1)-th\) 值执行除法来计算 \(\sigma\) 的值\((\sigma=\frac{M_{(1,i)}}{M_{(1,\ i-1)}})\)；
-
如果这不是第一次迭代 \((i>0)\)，则计算 \(\ \ \sigma\) 和 \(\sigma_{old}\) 之间的差值，即 \(∆=σ-σ_old\)；
-
如果 delta 值大于精度误差值 \(\varepsilon=10e-6\)\((e.g. ∆>ε)\)，则返回并继续步骤 2，否则终止计算并转到步骤 8；
-
完成上述计算后，我们将得到给定矩阵 \(A_{(mxm)}^\ast : \sigma_{max}\ \gets\sigma\) 的最大特征值；
下面的例子说明了以下计算
假设我们给定一个分解矩阵 \(A_{(mxm)}^\ast=A^TA\) 和单位向量 \(R=\left\{1,1,1\right\}\)
让我们通过执行以下计算来找到以下矩阵和向量的乘积
由于我们正在执行第一次迭代，且 \(i\ =\ 1\)，因此我们不执行其他计算，只是将获得的结果向量追加到矩阵 \(M\) 中。
之后，我们继续进行下一次迭代 \(i\ =\ 2\)。在这次迭代中，我们同样执行以下计算
同样，我们将中间结果附加到矩阵 \(M\) 中，并且，由于这不是第一次迭代，所以执行当前的特征值 \(\sigma\) 计算：\(\sigma=\frac{M_{(1,i)}}{M_{(1,\ i-1)}}=\frac{274}{17}=16.11\)。然后，我们计算 \(∆=σ-σ_old=16.11-0.00=16.11\) 的值。显然，\(∆\) 的值远大于精度误差 \(\varepsilon\) 的值。正因为如此，我们继续进行下一次迭代，在此期间我们将获得以下结果
之后，我们再次计算 \(∆\) 的值，并继续进行下一个计算阶段，直到 \(∆\) 的值实际小于给定的精度误差值 \(\varepsilon\)。完成上述计算后，我们将得到以下结果矩阵 \(M\) 和分解矩阵 \(A^\ast\) 的第一个最大特征值 \(\sigma_{max}\)
完成上述计算后，我们将得到以下结果矩阵 \(M\) 和分解矩阵 \(A^\ast\) 的第一个最大特征值 \(\sigma_{max}\)
| 17 | 274 | 4293 | .. | .. | .. | .. | .. | 3.9148848553274702e+28 | 6.0410605508676628e+29 |
| 16 | 241 | 3686 | .. | .. | .. | .. | .. | 3.3162556201511577e+28 | 5.1173155134159651e+29 |
| 12 | 174 | 2623 | .. | .. | .. | .. | .. | 2.3402394273230705e+28 | 3.6112244948026083e+29 |
在本例中，为了找到给定分解矩阵的第一个最大特征值，我们执行了 \(t = 25\) 次迭代，最终得到以下结果
步骤 2：计算矩阵 \(A^\ast\) 的特征向量
既然我们已经成功计算了第一个最大特征值 \(\sigma_{max}\)，现在让我们来找到与此特征值对应的特定特征向量。为此，我们必须将此特征值从给定矩阵 \(A^\ast=A^TA\) 的每个对角元素中减去，然后通过对给定矩阵执行 Jordan-Gaussian 变换以得到行阶梯形来找到线性方程组的非平凡解
此刻，让我们讨论如何使用 Jordan-Gaussian 变换来解决以下线性方程组，这是一个可以表述为以下计算算法的算法：
-
设 \(A_{(mxm)}^\ast={(A}^TA-\sigma I)=0\) – 是我们要解的线性方程组，\(i=\bar{1,m}\) – 是一个循环计数器变量；
-
将给定矩阵的 \(i-th\) 行作为主导行，并获得基准元素 \(\alpha=a_{i,i}\) 的值；
-
检查给定基准元素 \(\alpha\) 是否等于 0 或 1。如果不是，继续步骤 4，否则转到步骤 5；
-
将 \(i-th\) 主导行中的每个元素除以基准元素 \(\alpha\) 的值；
-
检查 \(\alpha\) 的值是否等于 0。如果是，交换 \(i-th\) 和 \((i+1)-th\) 行，否则转到下一步；
-
对于每个行 \(r=\bar{1,m}\)，除了主导行 \(i\)，通过重新计算每行的值来执行更新。为此，将每行的基准元素值作为 \(\gamma=a_{r,i}\)，然后使用以下公式更新这些值 \(k=\bar{1,m}\)： \(a_{r,k}=a_{r,k}-a_{i,k}\ast\ \gamma\)；
-
检查下一个基元素 \(\alpha_{i,i}\) 是否等于 0 或 \(i\ =\ m-1\)（例如，这是矩阵 \(A^\ast\) 的最后一列）。如果是，则跳到步骤 8，否则将 \(i\ =\ i\ +\ 1\) 的值递增并返回步骤 2 以继续下一次迭代。
-
从给定矩阵的 \(i\ +\ 1\ – th\) 列中提取每个值并将其赋给向量 \(X\)，该向量是给定线性方程组的非平凡解；
-
将向量 \(X\) 的每个分量值取反；
-
将向量 \(X\) 的最后一个分量赋值为 1，以获得与特征值 \(\sigma\) 完全对应的特征向量；
以下是执行此类计算的示例
假设我们已经给定一个在步骤 2 开始时计算的矩阵 \(A^\ast\)
让我们使用上面列出的算法来计算这个线性方程组的非平凡解。为此，我们取第一行作为主导行，并获得其基准元素 \(\alpha=a_{1,1}=-5.4310\)。此时，我们所要做的就是检查该值是否不等于 0 或 1，如果是，则将第一行中的每个元素除以 \(\alpha\) 的值
之后，我们必须通过执行以下计算来更新除主导行之外的每一行。对于每个特定行，我们必须获取基准元素 \(\gamma=a_{r,i}\) 的值。对于第二行，此值等于 \(\gamma=5\)。让我们将 \(\gamma\) 的值和主导第一行中的每个值应用于矩阵 \(A^\ast\) 的其他行
由于我们已经得到了第一主导行的值并更新了除第一行之外的每行的值，现在我们继续用第二主导行进行计算。为此，我们再次需要得到基本元素 \(\alpha=a_{2,2}=-4.8278\) 的值，并执行与上面已经完成的类似计算
与上一个阶段类似，我们还需要更新所有其他行中的值。为此，我们必须获取 \(\gamma=-0.9206\) 的值并执行以下计算
在计算了第二主导行并更新了其余行之后，我们可能会注意到矩阵 \(A^\ast\) 的所有对角元素都等于 1，除了最后一个对角元素 \(a_{2,2}\) 等于 0。这实际上意味着列 \(i=2\) 的所有元素都是这个方程组的非平凡解。此时，我们必须提取列 \(i=2\) 的所有元素并将它们放入向量 \(X\) 中，同时将每个元素取反。此外，我们必须将列中的最后一个元素（即向量 \(X\)）赋值为 1
显然，以下向量 X 是一个恰好对应于给定特征值 \(\sigma=15.4310\) 的特征向量。
步骤 3：寻找矩阵 \(A^\ast\) 的等价矩阵 \(B\)
由于我们已经计算出分解矩阵 \(A^\ast\) 的第一个最大特征值 \(\sigma_1\) 和特征向量 \(X_1\)，现在我们来寻找该矩阵的第二个特征值和特征向量。这通常通过计算一个尺寸减小的矩阵 \(B\) 来完成，该矩阵等价于矩阵 \(A^\ast\)，然后执行与步骤 1 和 2 中已经完成的相同计算，以找到给定矩阵 \(A^\ast\) 的第二对特征值和特征向量。为了找到等价矩阵 \(B\)，我们通常使用以下公式
其中 \(H\) – 矩阵 \(A^\ast\) 的厄米共轭，\(H^{-1}\) - 矩阵 \(A^\ast\) 的厄米共轭逆，\(B\) – 矩阵 \(A^\ast\) 的等价矩阵，计算其第二个特征值和特征向量；
为了找到等价矩阵 \(B\)，我们必须首先找到 Hermitian 共轭矩阵 \(H\) 及其逆 \(H^{-1}\)。 Hermitian 共轭矩阵 \(H\) 是对角单位矩阵的特例，可以表示如下：
Hermitian 矩阵 \(H\) 第一列中的每个元素都赋值为特征向量 \(x_i\) 的每个分量除以第一个分量 \(x_1\) 的分数形式，除了第一个元素是 \({1/x}_1\) 的除法。
此外，为了计算矩阵 \(H_{(mxm)}^{-1}\)，我们实际上不需要找到 Hermitian 矩阵 \(H\) 的逆。 Hermitian 共轭逆矩阵 \(H_{(mxm)}^{-1}\) 可以表示如下：
这个矩阵的第一列与Hermitian矩阵\(H\)非常相似，唯一的区别是我们将向量\(X\)的分量沿着第一列放置，除了第一个值外，其余值都取反，并且不将其除以第一个分量\(x_1\)。
以下示例演示了如何查找厄米特矩阵 \(H\) 及其逆 \(H^{-1}\)
既然我们已经找到了厄米特矩阵及其逆，现在让我们使用上面列出的公式计算等价矩阵 \(B\)
在执行了上述计算后，我们来寻找等价矩阵 \(B\)
最后，既然我们已经找到了等价矩阵 \(B\)，让我们通过执行与步骤 1 和 2 中已经完成的相同计算，来计算矩阵 \(A^\ast\gets B\) 的第二个特征值。
注意：在步骤2中计算出的每个特征向量实际上*不是*我们输入分解矩阵\(A^TA\)的特征向量。我们故意计算这些特征向量，以便能够通过Hermitian矩阵\(H\)及其逆矩阵\(H^{-1}\)的计算来找到等价矩阵\(B\)。对于每个等价矩阵\(B\)，在步骤1和2的每次迭代中，我们递归地旨在计算一个最大特征值\(\sigma_{max}\)，它是分解矩阵\(A^TA\)的下一个特征值。分解矩阵\(A^TA\)的特征向量在下一步骤4中单独计算，因为\(A^TA\)矩阵的所有特征值都已找到。
在步骤 1 和 2 中执行计算时，我们实际上正在寻找分解矩阵 \(A^TA\) 的每个特征值。我们执行以下计算，直到找到所有特定的特征值，然后继续执行步骤 4，以找到分解矩阵 \(A^TA\) 的所有特征向量。
这就是为什么，此时我们所要做的就是回到步骤 1，并继续计算分解矩阵 \(A^TA\) 的下一个特征值。
步骤 4：计算矩阵 \(A^TA\) 的特征向量
既然我们已经计算了矩阵 \(A^TA\) 的所有特征值 \((\sigma_1=15.4310,\ \sigma_2=5.5573,\ \sigma_3=0.0116)\)，现在我们可以使用步骤 2 中执行的 Jordan-Gaussian 变换，为分解矩阵 \(A^TA\) 的每个特征值找到一个特征向量。执行这些计算的结果，我们将得到以下特征向量
Using the Code
这是计算 SVD 的主函数
template<typename arg="long" double="">
void svd(std::vector<std::vector<arg>> matrix, std::vector<std::vector<arg>>& s,
std::vector<std::vector<arg>>& u, std::vector<std::vector<arg>>& v)
{
std::vector<std::vector<arg>> matrix_t;
matrix_transpose(matrix, matrix_t);
std::vector<std::vector<arg>> matrix_product1;
matrix_by_matrix(matrix, matrix_t, matrix_product1);
std::vector<std::vector<arg>> matrix_product2;
matrix_by_matrix(matrix_t, matrix, matrix_product2);
std::vector<std::vector<arg>> u_1;
std::vector<std::vector<arg>> v_1;
std::vector<arg> eigenvalues;
compute_evd(matrix_product2, eigenvalues, v_1, 0);
matrix_transpose(v_1, v);
s.resize(matrix.size());
for (std::uint32_t index = 0; index < eigenvalues.size(); index++)
{
s[index].resize(eigenvalues.size());
s[index][index] = eigenvalues[index];
}
std::vector<std::vector<arg>> s_inverse;
get_inverse_diagonal_matrix(s, s_inverse);
std::vector<std::vector<arg>> av_matrix;
matrix_by_matrix(matrix, v, av_matrix);
matrix_by_matrix(av_matrix, s_inverse, u);
}
此函数执行给定分解矩阵的特征值和特征向量的计算
template<typename arg="long" double="">
void compute_evd(std::vector<std::vector<arg>> matrix,
std::vector<arg>& eigenvalues, std::vector<std::vector<arg>>& eigenvectors,
std::size_t eig_count)
{
std::size_t m_size = matrix.size();
std::vector<arg> vec; vec.resize(m_size);
std::fill_n(vec.begin(), m_size, 1);
static std::vector<std::vector<arg>> matrix_i;
if (eigenvalues.size() == 0 && eigenvectors.size() == 0)
{
eigenvalues.resize(m_size);
eigenvectors.resize(eigenvalues.size());
matrix_i = matrix;
}
std::vector<std::vector<arg>> m; m.resize(m_size);
for (std::uint32_t row = 0; row < m_size; row++)
m[row].resize(100);
Arg lambda_old = 0;
std::uint32_t index = 0; bool is_eval = false;
while (is_eval == false)
{
for (std::uint32_t row = 0; row < m_size && (index % 100) == 0; row++)
m[row].resize(m[row].size() + 100);
for (std::uint32_t row = 0; row < m_size; row++)
{
m[row][index] = 0;
for (std::uint32_t col = 0; col < m_size; col++)
m[row][index] += matrix[row][col] * vec[col];
}
for (std::uint32_t col = 0; col < m_size; col++)
vec[col] = m[col][index];
if (index > 0)
{
Arg lambda = (m[0][index - 1] != 0) ? \
(m[0][index] / m[0][index - 1]) : m[0][index];
is_eval = (std::fabs(lambda - lambda_old) < 0.0000000001) ? true : false;
lambda = (std::fabs(lambda) >= 10e-6) ? lambda : 0;
eigenvalues[eig_count] = lambda; lambda_old = lambda;
}
index++;
}
std::vector<std::vector<arg>> matrix_new;
if (m_size > 1)
{
std::vector<std::vector<arg>> matrix_target;
matrix_target.resize(m_size);
for (std::uint32_t row = 0; row < m_size; row++)
{
matrix_target[row].resize(m_size);
for (std::uint32_t col = 0; col < m_size; col++)
matrix_target[row][col] = (row == col) ? \
(matrix[row][col] - eigenvalues[eig_count]) : matrix[row][col];
}
std::vector<arg> eigenvector;
jordan_gaussian_transform(matrix_target, eigenvector);
std::vector<std::vector<arg>> hermitian_matrix;
get_hermitian_matrix(eigenvector, hermitian_matrix);
std::vector<std::vector<arg>> ha_matrix_product;
matrix_by_matrix(hermitian_matrix, matrix, ha_matrix_product);
std::vector<std::vector<arg>> inverse_hermitian_matrix;
get_hermitian_matrix_inverse(eigenvector, inverse_hermitian_matrix);
std::vector<std::vector<arg>> iha_matrix_product;
matrix_by_matrix(ha_matrix_product,
inverse_hermitian_matrix, iha_matrix_product);
get_reduced_matrix(iha_matrix_product, matrix_new, m_size - 1);
}
if (m_size <= 1)
{
for (std::uint32_t index = 0; index < eigenvalues.size(); index++)
{
Arg lambda = eigenvalues[index];
std::vector<std::vector<arg>> matrix_target;
matrix_target.resize(matrix_i.size());
for (std::uint32_t row = 0; row < matrix_i.size(); row++)
{
matrix_target[row].resize(matrix_i.size());
for (std::uint32_t col = 0; col < matrix_i.size(); col++)
matrix_target[row][col] = (row == col) ? \
(matrix_i[row][col] - lambda) : matrix_i[row][col];
}
eigenvectors.resize(matrix_i.size());
jordan_gaussian_transform(matrix_target, eigenvectors[index]);
Arg eigsum_sq = 0;
for (std::uint32_t v = 0; v < eigenvectors[index].size(); v++)
eigsum_sq += std::pow(eigenvectors[index][v], 2.0);
for (std::uint32_t v = 0; v < eigenvectors[index].size(); v++)
eigenvectors[index][v] /= sqrt(eigsum_sq);
eigenvalues[index] = std::sqrt(eigenvalues[index]);
}
return;
}
compute_evd(matrix_new, eigenvalues, eigenvectors, eig_count + 1);
return;
}
这两个函数允许我们找到厄米特矩阵及其逆
template<typename arg="long" double="">
void get_hermitian_matrix(std::vector<arg> eigenvector,
std::vector<std::vector<arg>>& h_matrix)
{
h_matrix.resize(eigenvector.size());
for (std::uint32_t row = 0; row < eigenvector.size(); row++)
h_matrix[row].resize(eigenvector.size());
h_matrix[0][0] = 1 / eigenvector[0];
for (std::uint32_t row = 1; row < eigenvector.size(); row++)
h_matrix[row][0] = -eigenvector[row] / eigenvector[0];
for (std::uint32_t row = 1; row < eigenvector.size(); row++)
h_matrix[row][row] = 1;
}
template<typename arg="long" double="">
void get_hermitian_matrix_inverse(std::vector<arg> eigenvector,
std::vector<std::vector<arg>>& ih_matrix)
{
ih_matrix.resize(eigenvector.size());
for (std::uint32_t row = 0; row < eigenvector.size(); row++)
ih_matrix[row].resize(eigenvector.size());
ih_matrix[0][0] = eigenvector[0];
for (std::uint32_t row = 1; row < eigenvector.size(); row++)
ih_matrix[row][0] = -eigenvector[row];
for (std::uint32_t row = 1; row < eigenvector.size(); row++)
ih_matrix[row][row] = 1;
}
以下函数执行给定分解矩阵的 Jordan-高斯变换
template<typename arg="long" double="">
void jordan_gaussian_transform(
std::vector<std::vector<arg>> matrix, std::vector<arg>& eigenvector)
{
const Arg eps = 0.000001; bool eigenv_found = false;
for (std::uint32_t s = 0; s < matrix.size() - 1 && !eigenv_found; s++)
{
std::uint32_t col = s; Arg alpha = matrix[s][s];
while (col < matrix[s].size() && alpha != 0 && alpha != 1)
matrix[s][col++] /= alpha;
for (std::uint32_t col = s; col < matrix[s].size() && !alpha; col++)
std::swap(matrix[s][col], matrix[s + 1][col]);
for (std::uint32_t row = 0; row < matrix.size(); row++)
{
Arg gamma = matrix[row][s];
for (std::uint32_t col = s; col < matrix[row].size() && row != s; col++)
matrix[row][col] = matrix[row][col] - matrix[s][col] * gamma;
}
std::uint32_t row = 0;
while (row < matrix.size() &&
(s == matrix.size() - 1 || std::fabs(matrix[s + 1][s + 1]) < eps))
eigenvector.push_back(-matrix[row++][s + 1]);
if (eigenvector.size() == matrix.size())
{
eigenv_found = true; eigenvector[s + 1] = 1;
for (std::uint32_t index = s + 1; index < eigenvector.size(); index++)
eigenvector[index] = (std::fabs(eigenvector[index]) >= eps) ?
eigenvector[index] : 0;
}
}
}
以下函数用于寻找逆对角矩阵 S
template<typename arg="long" double="">
void get_inverse_diagonal_matrix(std::vector<std::vector<arg>> matrix,
std::vector<std::vector<arg>>& inv_matrix)
{
inv_matrix.resize(matrix.size());
for (std::uint32_t index = 0; index < matrix.size(); index++)
{
inv_matrix[index].resize(matrix[index].size());
inv_matrix[index][index] = 1.0 / matrix[index][index];
}
}
此函数允许我们计算等价矩阵 B
template<typename arg="long" double="">
void get_reduced_matrix(std::vector<std::vector<arg>> matrix,
std::vector<std::vector<arg>>& r_matrix, std::size_t new_size)
{
r_matrix.resize(new_size);
std::size_t index_d = matrix.size() - new_size;
std::uint32_t row = index_d, row_n = 0;
while (row < matrix.size())
{
r_matrix[row_n].resize(new_size);
std::uint32_t col = index_d, col_n = 0;
while (col < matrix.size())
r_matrix[row_n][col_n++] = matrix[row][col++];
row++; row_n++;
}
}
以下简单函数执行两个矩阵的乘法
template<typename arg="long" double="">
void matrix_by_matrix(std::vector<std::vector<arg>> matrix1,
std::vector<std::vector<arg>>& matrix2, std::vector<std::vector<arg>>& matrix3)
{
matrix3.resize(matrix1.size());
for (std::uint32_t row = 0; row < matrix1.size(); row++)
{
matrix3[row].resize(matrix1[row].size());
for (std::uint32_t col = 0; col < matrix1[row].size(); col++)
{
matrix3[row][col] = 0.00;
for (std::uint32_t k = 0; k < matrix1[row].size(); k++)
matrix3[row][col] += matrix1[row][k] * matrix2[k][col];
}
}
}
以下简单函数允许我们找到给定矩阵的转置
template<typename arg="long" double="">
void matrix_transpose(std::vector<std::vector<arg>> matrix1,
std::vector<std::vector<arg>>& matrix2)
{
matrix2.resize(matrix1.size());
for (std::uint32_t row = 0; row < matrix1.size(); row++)
{
matrix2[row].resize(matrix1[row].size());
for (std::uint32_t col = 0; col < matrix1[row].size(); col++)
matrix2[row][col] = matrix1[col][row];
}
}
</std::vector<arg></std::vector<arg></typename>
这两个函数允许我们生成并打印矩阵
template<typename arg="long" double="">
void generate_matrix(std::vector<std::vector<long double="">>& \
matrix, std::size_t rows, std::size_t cols)
{
std::srand((unsigned int)std::time(nullptr)); matrix.resize(rows);
for (std::size_t row = 0; row < matrix.size(); row++)
{
matrix[row].resize(cols);
for (std::size_t col = 0; col < matrix[row].size(); col++)
matrix[row][col] = std::rand() % 20 - std::rand() % 20;
}
}
template<typename arg="long" double="">
void print_matrix(std::vector<std::vector<long double="">> matrix)
{
for (std::size_t row = 0; row < matrix.size(); row++)
{
for (std::size_t col = 0; col < matrix[row].size(); col++)
std::cout << std::setprecision(5) << std::fixed << matrix[row][col] << " ";
std::cout << "\n";
}
std::cout << "\n";
}
以下主函数演示了完整的SVD计算
int main(int argc, char* argv[])
{
std::size_t matrix_size = 0;
std::vector<std::vector<long double="">> u, v;
std::vector<std::vector<long double="">> matrix, s;
std::cout << "Singular Value Decomposition (SVD):\n\n";
std::cout << "Enter size of matrix N = "; std::cin >> matrix_size;
generate_matrix(matrix, matrix_size, matrix_size);
std::cout << "\nA = \n"; print_matrix(matrix);
svd(matrix, s, u, v);
std::cout << "\nS = \n"; print_matrix(s);
std::cout << "\nU = \n"; print_matrix(u);
std::cout << "\nV = \n"; print_matrix(v);
std::cin.get();
std::cin.get();
return 0;
}
关注点
如今，奇异值分解主要用作潜在语义分析（LSA）的一部分。完整的 SVD 作为一种方法，允许我们将给定的实矩阵或酉矩阵表示在不同的正交空间（即坐标系）中。完整的 SVD 最适用于自然语言处理和分布语义，因为它允许我们检测给定频率事件矩阵的主成分，忽略可能的“噪声”。具体来说，如果我们有一个事件矩阵，其中每个元素是某个文档中每个短语出现的频率，那么 SVD 方法就非常有用。在这种情况下，计算完整的 SVD 允许我们将每个短语和包含该短语的文档分组，然后执行简单的聚类以找到包含某个短语的所有文档。
历史
- 2018年11月30日——本文第一版发布