问关于低秩矩阵的计算

EN

如果您知道矩阵的秩为3，那么恰好3次Householder变换就足以将nxm矩阵缩减为3xm矩阵。现在可以很容易地将其转换为特征值问题。计算特征多项式。例如，考虑这个矩阵(我将使用MATLAB来完成此工作)：

>> A = rand(10,3);
>> B = A*rand(3,6);

很明显，B的排名是3，但是如果你不相信我，排名肯定了这个说法。

>> rank(B)
ans =
     3

>> size(B)
ans =
    10     6

所以B是一个10x6的矩阵，秩为3。奇异值分解也证实了这一事实。

>> svd(B)
ans =
          6.95958965358531
          1.05904552889497
         0.505730235622534
      7.37626877572817e-16
      2.36322066654691e-16
      1.06396598411356e-16

我觉得太懒了，不想写“户主转换”。(我已经准备了一些代码，但是...)因此，我将使用QR来帮助我。

>> [Q,R] = qr(B,0);
>> C = Q(:,1:3)'*B
C =
   -2.0815   -1.7098   -3.7897   -1.6186   -3.6038   -3.0809
    0.0000    0.91329   0.78347   0.44597  -0.072369   0.54196
    0.0000    0.0000   -0.2285   -0.43721  -0.85949  -0.41072

这里的乘法显示了我们在三次Householder变换后会看到的情况。正如我们所期望的，它是上三角形的。接下来，计算特征多项式。我将在这里象征性地使用我自己的工具，但计算只是一点代数。

>> sympoly lambda
>> det(C*C' - lambda*eye(3))
ans =
    13.8942 + 66.9996*lambda + 49.8132*lambda^2 + lambda^3

>> P = det(C*C' - lambda*eye(3))
P =
    13.8942 - 66.9996*lambda + 49.8132*lambda^2 - lambda^3

P的根是什么，那么C*C‘的特征值是什么？

>> r = roots(P)
r =
       48.436
       1.1216
      0.25576

我们知道特征值必须是这里奇异值的平方，所以这里的一个很好的测试是看看我们是否恢复了svd找到的奇异值。因此，再次扩展显示格式，我们可以看到它做得非常好。

>> sqrt(roots(P))
ans =
          6.95958965358531
          1.05904552889497
         0.505730235622533

数学可以是有趣的，当它工作。我们能用这些信息做什么呢？如果我知道一个特定的特征值，我就可以计算相应的特征向量。本质上，我们需要求解线性3x3齐次方程组。

(C*C' - eye(3)*r(3)) * X = 0

再说一次，我懒得在这里找到解决方案，而不需要实际编写任何代码。高斯消去法可以做到这一点。

>> V = null((C*C' - eye(3)*r(3)))
V =
        -0.171504758161731
        -0.389921448437349
         0.904736084157367

所以我们有V，C*C‘的特征向量。我们可以通过以下对svd的调用来说服自己。

>> svd(C - V*(V'*C))
ans =
           6.9595896535853
          1.05904552889497
      2.88098729108798e-16

减去C在V方向上的分量，我们得到一个秩为2的矩阵。

类似地，我们可以将V转换到原始问题空间，并使用它来转换矩阵B，我们的原始矩阵，通过B的秩1更新。

>> U = Q(:,1:3)*V;
>> D = B - U*(U'*B);

D的排名是什么？

>> svd(D)
ans =
          6.95958965358531
          1.05904552889497
      2.62044567948618e-15
      3.18063391331806e-16
      2.16520878207897e-16
      1.56387805987859e-16

>> rank(D)
ans =
     2

正如你所看到的，尽管我做了很多数学计算，多次调用svd，QR，rank等，但最终，实际的计算是相对微不足道的。我只需要...

1. 3户主转换。(将它们存储起来以供以后使用。请注意，霍斯霍尔德变换只是对特征向量( matrix.)
2. Computate polynomial.
3. Compute )的秩1更新，对立方polynomial.
4. Recover (相应的特征向量)使用您最喜欢的方法。这里高斯消去法就足够了。
5. 使用我们最初做过的householder变换将特征向量返回到原始空间。
6. 对矩阵进行排名一的更新。

所有这些计算步骤对于任何大小的矩阵都是快速和有效的，只要我们知道实际的秩是3。甚至不值得在这个主题上发表论文。

编辑：

由于这个问题已经被修改，使得矩阵本身只有3x3的大小，所以计算甚至更简单。但是，我将保留这篇文章的第一部分，因为它描述了一个完全有效的解决方案，适用于任意大小的秩为3的通用矩阵。

如果你的目标是消除3x3矩阵的最后一个奇异值，那么3x3矩阵上的svd似乎是相当有效的。在用间接方法生成最后一个奇异值时，也会有一些精度损失。例如，在这里比较由svd计算的最小奇异值，然后使用特征值。所以你可能会在这里看到一些小错误，因为形成A'*A会导致一些精度。这种损失的程度可能取决于A的条件数。

>> A = rand(3);
>> sqrt(eig(A'*A))
ans =
        0.0138948003095069
         0.080275195586312
          1.50987693453097
>> svd(A)
ans =
          1.50987693453097
        0.0802751955863124
        0.0138948003095054

但是，如果你真的想自己做这项工作，你可以自己做。

1. 直接从3x3矩阵A'*A.
2. Compute该三次多项式的根计算特征多项式。在A上选择最小的root.
3. Generate对应的eigenvector.
4. Do进行秩1更新，去掉A中位于该eigenvector.

方向的那部分

与简单地调用svd，然后进行一次排名更新相比，这是否更简单或计算效率更低？我根本不确定是否值得在3x3上付出努力。一个3x3的svd计算起来真的非常快。

你可能已经考虑到了，但是如果A是正常的，SVD可以通过特征值分解来计算。这相当于求解特征多项式，并且因为它是秩3矩阵的3次多项式，所以根是由众所周知的立方公式找到的。

它看起来通常SVD必须简化为求解秩为3的矩阵的立方，但我不记得读过任何关于这方面的内容。一个谷歌快速转向的this piece of code，声称以这种方式解决了排名3的奇异值分解。不幸的是，没有随附的文件。如果使用此代码，使用非正定矩阵测试它应该是一个很好的测试用例。

编辑:在第二次阅读时，作者似乎也使用了特征分解。在非PD矩阵上可能不是很好，但我希望在这里证明我是错的。

您可以使用幂迭代来查找与最大奇异值对应的奇异向量。一旦你有了它，你就可以使用幂迭代的修改版本来找到第二大奇异值的向量；在每次迭代之后，你需要减去部分解向量到最大奇异值向量的投影。

这是一种寻找所有奇异向量的糟糕方法，但对于前两个应该可以很好地工作。