CODEWITHSUNDEEP
pythonnumpylinear-algebra
Seysen
Seysen

Reputation: 134

Performing Householder Reflection of a vector for QR Decomposition

This question was asked before me on here.

However, the solution there was not satisfactory for me, I am still stuck at 33% mismatch, so I felt the need to re-open this topic (And also the author of that thread didn't add an appropriate answer after solving the issue for themselves).

The code that I have written is here:

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

Here is the test case that this code is supposed to pass:

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

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

The first assertion is passed successfully, however, the second one gives me a 33% mismatch:

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.      ])

I have been trying everything for like 5 hours now, and I feel like I'm wasting too much time on this. Any help to make this code pass the test would be much appreciated by me.

Upvotes: 0

Views: 1578

Answers (1)

hfhc2
hfhc2

Reputation: 4391

Well, it looks correct to me.

The problem seem to be the parameters of the assert_allclose function. Specifically, it reports whether or not 

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

for each pair of entries a and b. According to the docs, the absolute tolerance is 1e-8 for the ordinary allclose function. However, the assert_allclose parameter of atol is 0 by default.

Since your target b is zero, any value != 0 is not close with respect to this function, even though the two values are certainly reasonably close.

I recommend setting atol to 1e-8, i.e. 

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

I am not quite sure why the numpy people chose different parameters for the ordinary allclose and assert_allclose though...

Upvotes: 1

Related Questions