【问题标题】:Calculating the null space of a matrix计算矩阵的零空间
【发布时间】:2011-02-28 21:27:18
【问题描述】:

我正在尝试求解一组 Ax = 0 形式的方程。A 是已知的 6x6 矩阵,我使用 SVD 编写了以下代码来获得在一定程度上起作用的向量 x。答案大致正确但不足以对我有用,如何提高计算的精度?将 eps 降低到 1.e-4 以下会导致函数失败。

from numpy.linalg import *
from numpy import *

A = matrix([[0.624010149127497 ,0.020915658603923 ,0.838082638087629 ,62.0778180312547 ,-0.336 ,0],
[0.669649399820597 ,0.344105317421833 ,0.0543868015800246 ,49.0194290212841 ,-0.267 ,0],
[0.473153758252885 ,0.366893577716959 ,0.924972565581684 ,186.071352614705 ,-1 ,0],
[0.0759305208803158 ,0.356365401030535 ,0.126682113674883 ,175.292109352674 ,0 ,-5.201],
[0.91160934274653 ,0.32447818779582 ,0.741382053883291 ,0.11536775372698 ,0 ,-0.034],
[0.480860406786873 ,0.903499596111067 ,0.542581424762866 ,32.782593418975 ,0 ,-1]])

def null(A, eps=1e-3):
  u,s,vh = svd(A,full_matrices=1,compute_uv=1)
  null_space = compress(s <= eps, vh, axis=0)
  return null_space.T

NS = null(A)
print "Null space equals ",NS,"\n"
print dot(A,NS)

【问题讨论】:

    标签: python math linear-algebra svd least-squares


    【解决方案1】:

    注意:python 与 matlab-syntax(?) 中的 SVD 可能存在混淆: 在 python 中,numpy.linalg.svd(A) 返回矩阵 u,s,v 使得 u*s*v = A (严格来说:dot(u, dot(diag(s), v) = A,因为 s 是向量而不是 numpy 中的二维矩阵)。

    从这个意义上说,最上面的答案是正确的,通常你写 u*s*vh = A 并返回 vh,这个答案讨论的是 v 而不是 vh。

    长话短说:如果你有矩阵 u,s,v 使得 u*s*v = A,那么 v 的最后 不是 v 的最后一列,描述零空间。

    编辑:[对于像我这样的人:]最后一行是一个向量 v0,使得 A*v0 = 0(如果相应的奇异值为 0)

    【讨论】:

      【解决方案2】:

      A 是满级 --- 所以x0

      因为看起来您需要一个最小二乘解决方案,即min ||A*x|| s.t. ||x|| = 1,所以执行 SVD 使得[U S V] = svd(A)V 的最后一列(假设这些列按奇异值递减的顺序排序)是x

      即,

      U =
      
           -0.23024     -0.23241      0.28225     -0.59968     -0.04403     -0.67213
            -0.1818     -0.16426      0.18132      0.39639      0.83929     -0.21343
           -0.69008     -0.59685     -0.18202      0.10908     -0.20664      0.28255
           -0.65033      0.73984    -0.066702     -0.12447     0.088364       0.0442
        -0.00045131    -0.043887      0.71552     -0.32745       0.1436      0.59855
           -0.12164      0.11611       0.5813      0.59046     -0.47173     -0.25029
      
      
      S =
      
             269.62            0            0            0            0            0
                  0       4.1038            0            0            0            0
                  0            0        1.656            0            0            0
                  0            0            0       0.6416            0            0
                  0            0            0            0      0.49215            0
                  0            0            0            0            0   0.00027528
      
      
      V =
      
          -0.002597     -0.11341      0.68728     -0.12654      0.70622    0.0050325
         -0.0024567     0.018021       0.4439      0.85217     -0.27644    0.0028357
         -0.0036713      -0.1539      0.55281      -0.4961      -0.6516   0.00013067
            -0.9999    -0.011204   -0.0068651    0.0013713    0.0014128    0.0052698
          0.0030264      0.17515      0.02341    -0.020917   -0.0054032      0.98402
           0.012996     -0.96557     -0.15623      0.10603     0.014754      0.17788
      

      所以,

      x =
      
          0.0050325
          0.0028357
         0.00013067
          0.0052698
            0.98402
            0.17788
      

      并且,||A*x|| = 0.00027528 与您之前针对 x 的解决方案相反,其中 ||A*x_old|| = 0.079442

      【讨论】:

      • x=0 是解决问题的方法,但也很无趣。通过不同方式得出的问题的真正解决方案是:[0.880057009282733,0.571293018023548,0.0664250041765576,1,186.758799941964,33.7579819749057]T
      • 你确定吗?我在 A*x 的结果中看到了一些非零元素 --- [-0.056356 -0.055643 -7.3896e-013 -0.0043278 0.004483 -2.1316e-014]
      • 当然,除非你不想要零空间,而是最小二乘解决方案,即min ||A*x|| s.t. ||x|| = 1
      • 我同意雅各布的观点。 A 有满秩。 eps 为 1e-4 时出现错误的原因是矩阵的最小奇异值是 2.75282332e-04。换句话说,您需要具有 0 的奇异值(在浮点精度范围内)才能具有包含零向量以外的向量的零空间。顺便说一句,Matlab 也将x 设为 0。
      • 更新了最小二乘解。
      猜你喜欢
      • 1970-01-01
      • 2018-02-25
      • 2023-04-02
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2016-01-29
      • 1970-01-01
      相关资源
      最近更新 更多