0

I am trying to implement the QR decomposition via householder reflectors. While attempting this on a very simple array, I am getting weird numbers. Anyone who can tell me, also, why using the @ vs * operator between vec and vec.T on the last line of the function definition gets major bonus points.

This has stumped two math/comp sci phds as of this morning.

    import numpy as np

    def householder(vec):
       vec[0] += np.sign(vec[0])*np.linalg.norm(vec)
       vec = vec/vec[0]
       gamma = 2/(np.linalg.norm(vec)**2)
       return np.identity(len(vec)) - gamma*(vec*vec.T)

    array = np.array([1, 3 ,4])
    Q = householder(array)
    print(Q@array)

Output: 

   array([-4.06557377, -7.06557377, -6.06557377])

Where it should be:

   array([5.09, 0, 0])
F.Dunbar
  • 37
  • 6

1 Answers1

0

* is elementwise multiplication, @ is matrix multiplication. Both have their uses, but for matrix calculations you most likely want the matrix product.

vec.T for an array returns the same array. A simple array only has one dimension, there is nothing to transpose. vec*vec.T just returns the elementwise squared array.

You might want to use vec=vec.reshape(-1,1) to get a proper column vector, a one-column matrix. Then vec*vec.T does "by accident" the correct thing. You might want to put the matrix multiplication operator there anyway.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51