Optimized projection of a matrix orthogonaly to a vector with Numpy

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.

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