假设我有一个秩为 3 的矩阵 3x3 A。 我想创建一个在 $ {l}_{2} $ / Frobenius 范数中最接近 A 的 Rank 2 矩阵。 我们称这个矩阵为 F。
由SVD很容易实现，即如果$A=U S {V}^{H}$由SVD分解$F=U \hat{S} {V}^{H}$。 其中 $ \hat{S} $ 与 $ S $ 相同，最后一个奇异值为零。
问题是，有没有一种计算量较小的方法来创建 F 但使用 SVD 分解？
谢谢。
【问题讨论】：
标签： performance algorithm math matrix
如果您知道矩阵的秩为 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 也证实了这一事实。
>> svd(B)
ans =
6.95958965358531
1.05904552889497
0.505730235622534
7.37626877572817e-16
2.36322066654691e-16
1.06396598411356e-16
我懒得写 Householder 转换。 （我有一些代码，但是......）所以我会使用 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
因此，通过在 V 方向上减去 C 的那个分量，我们得到一个秩为 2 的矩阵。
类似地，我们可以将 V 转换为原始问题空间，并使用它通过 B 的一级更新来转换矩阵 B，即我们的原始矩阵。
>> 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 等，但最终实际计算还是比较琐碎的。我只需要...
只要我们知道实际排名为 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
但是，如果你真的想自己做这项工作，你可以做到。
这是否比简单调用 svd 更简单或计算效率更低，然后进行一级更新？我完全不确定在 3x3 上付出努力是否值得。 3x3 svd 的计算速度确实非常快。
【讨论】：
您可以使用幂迭代来找到最大奇异值对应的奇异向量。一旦你有了它，你可以使用修改版本的幂迭代来找到第二大奇异值的向量；每次迭代后，您需要将部分解向量的投影减去最大奇异值的向量。
这是查找所有奇异向量的糟糕方法，但对于前两个向量应该可以正常工作。
【讨论】：
您可能已经考虑过，但是如果 A 是正态的，则可以通过特征值分解来计算 SVD。这相当于求解特征多项式，并且由于它是 3 阶矩阵的 3 次多项式，因此可以通过众所周知的三次公式找到根。
看起来一般情况下 SVD 必须可简化为求解 3 阶矩阵的三次方，但我不记得读过任何关于此的内容。一个快速的 google-around 转向 this piece of code，声称以这种方式解决了 3 级 SVD。不幸的是，没有随附的文件。如果使用此代码，使用非正定矩阵对其进行测试应该是一个很好的测试用例。
编辑：在二读时，似乎那里的作者也使用了特征分解。在非 PD 矩阵上可能不好，但我想在这里被证明是错误的。
【讨论】：