Generate a vector that is orthogonal to a set of other vectors in any dimension

Assume I have a set of vectors $ a_1, ..., a_d $ that are orthonormal to each other. Now, I want to find another vector $ a_{d+1} $ that is orthogonal to all the other vectors.

1. Is there an efficient algorithm to achieve this? I can only think of adding a random vector to the end, and then applying gram-schmidt.

2. Is there a python library which already achieves this?

Related. Can't speak to optimality, but here is a working solution. The good thing is that numpy.linalg does all of the heavy lifting, so this may be speedier and more robust than doing Gram-Schmidt by hand. Besides, this suggests that the complexity is not worse than Gram-Schmidt.

The idea:

1. Treat your input orthogonal vectors as columns of a matrix O.
2. Add another random column to O. Generically O will remain a full-rank matrix.
3. Choose b = [0, 0, ..., 0, 1] with len(b) = d + 1.
4. Solve a least-squares problem x O = b. Then, x is guaranteed to be non-zero and orthogonal to all original columns of O.

import numpy as np
from numpy.linalg import lstsq
from scipy.linalg import orth

# random matrix
M = np.random.rand(10, 5)

# get 5 orthogonal vectors in 10 dimensions in a matrix form
O = orth(M)


def find_orth(O):
    rand_vec = np.random.rand(O.shape[0], 1)
    A = np.hstack((O, rand_vec))
    b = np.zeros(O.shape[1] + 1)
    b[-1] = 1
    return lstsq(A.T, b)[0]


res = find_orth(O)

if all(np.abs(np.dot(res, col)) < 10e-9 for col in O.T):
    print("Success")
else:
    print("Failure")