当前位置：文江博客话题详情

低阶矩阵的计算

发布于 2025-01-03 08:30:52 字数 248 浏览 0 评论 0原文

假设我有一个 3x3 A 的矩阵 3。我想创建一个 Rank 2 矩阵，它在 $ {l}_{2} $ / Frobenius 范数中最接近 A。我们称这个矩阵为F。

通过SVD很容易实现，即如果$A=US{V}^{H}$通过SVD分解$F=U\hat{S}{V}^{H}$ 。其中 $ \hat{S} $ 与 $ S $ 相同，最后一个奇异值为零。

问题是，是否有一种计算强度较小但使用 SVD 分解来创建 F 的方法？

谢谢。

原文

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

╰つ倒转 2025-01-10 08:30:52

如果您知道矩阵的秩为 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

因此，通过减去 C 在 V 方向上的分量，我们得到一个 2 阶矩阵。

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

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

的秩是多少？

>> 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 次户主转变。（存储它们以供以后使用。请注意，Householder 变换只是对矩阵的一级更新。）
计算特征多项式。
使用您最喜欢的三次多项式方法计算最小根。
恢复相应的特征向量。这里高斯消除就足够了。
使用我们最初所做的 Housener 变换将该特征向量返回到原始空间。
对矩阵的一级更新。

只要我们知道实际的秩是 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

但是，如果您确实想自己做这项工作，您可以做。

直接根据 3x3 矩阵 A'*A 计算特征多项式。
计算该三次多项式的根。选择最小的根。
生成相应的特征向量。
对 A 进行一级更新，删除 A 位于该特征向量方向的部分。

与简单调用 svd，然后进行一级更新相比，这是否更简单或计算效率更低？我完全不确定在 3x3 上付出努力是否值得。 3x3 svd 的计算速度确实非常快。

If you KNOW the matrix to be rank 3, then exactly 3 Householder transformations will be sufficient to reduce an nxm matrix to a 3xm matrix. One can now easily convert this into an eigenvalue problem. Compute the characteristic polynomial. For example, consider this matrix (I'll use MATLAB to do the work):

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

Clearly, B will have rank 3, but if you don't trust me, rank confirms that claim.

>> rank(B)
ans =
     3

>> size(B)
ans =
    10     6

So B is a 10x6 matrix, with rank 3. The SVD confirms that fact also.

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

I'm feeling too lazy to write the Householder transformations. (I've got some code for that laying around, but...) So I'll use QR to help me out.

>> [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

The multiply here shows what we would have seen after three Householder transformations. See that it is upper triangular as we would expect. Next, compute the characteristic polynomial. I'll so it symbolically here, using my own tools, but the computation is just a bit of algebra.

>> 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

What are the roots of P, so the eigenvalues of C*C'?

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

We know that the eigenvalues must be the squares of the singular values here, so a good test here is to see if we recovered the singlular values that svd found. So expanding the display format again, we see it did so nicely.

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

Mathematics can be fun when it works. What can we do with this information? If I know a specific eigenvalue, I can compute the corresponding eigenvector. Essentially, we need to solve the linear 3x3 homogeneous system of equations

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

Again, I'll be lazy here to find the solution without actually writing any code. Gaussian elimination would do it though.

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

So we have V, an eigenvector of C*C'. We can convince ourselves that it is one by the following call to svd.

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

So by subtracting off that component of C in the direction of V, we get a rank 2 matrix.

Similarly, we can transform V into the original problem space, and use that to transform the matrix B, our original matrix, by a rank one update of B.

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

What is the rank of D?

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

>> rank(D)
ans =
     2

As you can see, even though I did a lot of mathematics, calling svd, QR, rank, etc., all multiple times, in the end, the actual computation was relatively trivial. I needed only...

3 householder transformations. (Store them for later use. Note that a Householder transformation is simply a rank one update to your matrix.)
Computate the characteristic polynomial.
Compute the smallest root using your favorite method for a cubic polynomial.
Recover the corresponding eigenvector. Gaussian elimination is sufficient here.
Back that eigenvector into the original space using the householder transformations we originally did.
A rank one update to the matrix.

All of these computational steps would be fast and efficient for ANY size matrix, as long as we know the actual rank is 3. Not even worth a paper on the subject.

EDIT:

Since the question has been revised such that the matrix itself is only of size 3x3, the computation is even more trivial. I'll leave the first part of the post as is however, since it describes a fully valid solution for a general matrix of any size that is of rank 3.

If your goal is to kill off the last singular value of a 3x3 matrix, the svd on a 3x3 seems like it would be pretty efficient to do. There will also be some loss of precision in generating the last singular value by indirect means. For example, compare here the smallest singular value computed by svd, then by using the eigenvalues. So you may see some small error here, since forming A'*A will cause some of precision. The extent of this loss will probably depend on the condition number of A.

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

However, if you really want to do the work yourself, you could do it.

Compute the characteristic polynomial directly from the 3x3 matrix A'*A.
Compute the roots of that cubic polynomial. Choose the smallest root.
Generate the corresponding eigenvector.
Do a rank one update on A, killing off the part of A that lies in the direction of that eigenvector.

Is this any simpler or less computationally efficient than a simple call to the svd, followed by a rank one update? I'm not at all sure it is worth the effort on a 3x3. A 3x3 svd really is extremely fast to compute.

回复收藏 0 原文