执行向量的 Householder Reflection 以进行 QR 分解

Question

我在 here 上问过这个问题。

然而，那里的解决方案并不令我满意，我仍然停留在 33% 不匹配，所以我觉得有必要重新打开这个主题（而且该线程的作者没有添加适当的自己解决问题后再回答）。

我写的代码在这里：

def householder(vec):
    vec = np.asarray(vec, dtype=float)
    if vec.ndim != 1:
        raise ValueError("vec.ndim = %s, expected 1" % vec.ndim)

    n = len(vec)
    I = np.eye(n)
    e1 = np.zeros_like(vec).astype(float)
    e1[0] = 1.0
    V1 = e1 * np.linalg.norm(vec)
    print("V1:", V1)
    u = vec
    u[0] = -(np.sum(np.square(u[1:]))) / (vec[0] + np.linalg.norm(vec))
    u = u / np.linalg.norm(u)
    H = I - 2 * (np.outer(u, u))
    return V1 , H

下面是这段代码应该通过的测试用例：

v = np.array([1, 2, 3])
v1, h = householder(v)

assert_allclose(np.dot(h, v1), v)
assert_allclose(np.dot(h, v), v1)

第一个断言成功通过，但是，第二个断言给了我 33% 的不匹配：

AssertionError: 
Not equal to tolerance rtol=1e-07, atol=0

Mismatch: 33.3%
Max absolute difference: 4.4408921e-16
Max relative difference: 1.18687834e-16
 x: array([3.741657e+00, 2.220446e-16, 0.000000e+00])
 y: array([3.741657, 0.      , 0.      ])

我已经尝试了大约 5 个小时的所有方法，我觉得我在这上面浪费了太多时间。任何帮助使此代码通过测试的帮助我都将不胜感激。

Answer 1

嗯，我觉得是对的。

问题似乎是 assert_allclose 函数的参数。具体来说，它报告是否

absolute(a - b) <= (atol + rtol * absolute(b))

每对条目 a 和 b。根据docs, the absolute tolerance is 1e-8 for the ordinary allclose function. However, the assert_allclose parameter of atol is 0 by default。

由于您的目标 b 为零，任何值 != 0 相对于此函数都不接近，即使这两个值肯定 合理地 接近。

我建议将 atol 设置为 1e-8，即

assert_allclose(np.dot(h, v), v1, atol=1e-8)

我不太清楚为什么 numpy 的人会为普通 allclose 和 assert_allclose 选择不同的参数...

执行向量的 Householder Reflection 以进行 QR 分解

Performing Householder Reflection of a vector for QR Decomposition

python

numpy

linear-algebra