0

I need to make all other columns of a matrix A orthogonal to one of its column j.

I use the following algorithm :

# Orthogonalize with selected column
for i in remaining_cols:
    A[:,i] = A[:,i] - A[:,j] * np.dot(A[:,i], A[:,j]) / np.sum(A[:,j]**2)

The idea comes from the QR decomposition with the Gram-Schmidt process.

But this code is not optimized and unstable because of the Gram-Schmidt process.

Does Numpy provide any method to compute the orthogonal projection of those vectors ?

With Householder Matrix

I heard that the Householder Reflectors are used in numpy.linalg.qr. This would allow me to compute an orthogonal matrix Q so that 

Q * A[:,j] = [0 ... 0 1 0 ... 0]
                      |
                 j_th coordinate

I would only have to ignore the line j and multiply back with Q.T.

Is there a method to obtain the Householder Matrix with Numpy ? I mean without coding the algorithm by hand.

Matthias Beaupère
  • 1,731
  • 2
  • 17
  • 44
  • Have read this [post](https://stackoverflow.com/questions/17836880/orthogonal-projection-with-numpy)? – Ben.T May 28 '19 at 13:36
  • @Ben.T I did come across it, but I thought it was unrelated. Indeed I don't deal with Least Square problem, but maybe there is a connection ? – Matthias Beaupère May 28 '19 at 13:43
  • 1
    @Ben.T to my knowledge, LS problem enables to find an hyperplane closest to your data in terms of the euclidean distance. In my case, I have a specific projection hyperplane, the one that has the `A[:,j]`as normal vector. – Matthias Beaupère May 28 '19 at 14:02

1 Answers1

1

IIUC, here could be a vectorized way:

np.random.seed(10)
B = np.random.rand(3,3)

col = 0
remaining_cols = [1,2]

#your method
A = B.copy()
for i in remaining_cols:
    A[:,i] = A[:,i] - A[:,col] * np.dot(A[:,i], A[:,col]) / np.sum(A[:,col]**2)
print (A)
[[ 0.77132064 -0.32778252  0.18786796]
 [ 0.74880388  0.16014712 -0.2079702 ]
 [ 0.19806286  0.67103261  0.05464156]]

# vectorize method
A = B.copy()
A[:,remaining_cols] -= (A[:,col][:,None] * np.sum(A[:,remaining_cols]* A[:,col][:,None], axis=0)
                                              / np.sum(A[:,col]**2))

print (A) #same result
[[ 0.77132064 -0.32778252  0.18786796]
 [ 0.74880388  0.16014712 -0.2079702 ]
 [ 0.19806286  0.67103261  0.05464156]]
Ben.T
  • 29,160
  • 6
  • 32
  • 54
  • Works pretty nicely thanks ! I haven't accepted to see if there can be other propositions, especially using some Numpy optimized method. – Matthias Beaupère May 29 '19 at 14:34
  • @matthiasbe good to know it works and no pb, I assume there is an implemented method, but broadcasting is usually good regarding performance as far I know :) – Ben.T May 29 '19 at 15:18