I have a couple of orthonormal vectors. I would like to extend this 2-dimensional basis to a larger one. What is the fastest way of doing this in Python with NumPy?

My thoughts were the following: Generate a random vector of the required size (new_dimension > 2), perform Gram-Schmidt by substracting scaled dot-products with the previous two. Repeat. I doubt that this is the quickest way though...

CodePudding user response：

You didn't specify the dimension of your space. If it is 3, then you can simply use the cross product of your two vectors. If it is not, then see below.

Example in 3-D

# 1. setup: an orthonormal basis of two vectors a, b
np.random.seed(0)
a, b = np.random.uniform(size=(2,3))
a /= np.linalg.norm(a)
b -= a.dot(b)*a
b /= np.linalg.norm(b)

# 2. check:
>>> np.allclose([1,1,0,0], [a.dot(a), b.dot(b), a.dot(b), b.dot(a)])
True

Then, making a new vector:

# 3. solve
c = np.cross(a, b)

# 4. checks
>>> np.allclose([1,0,0], [c.dot(c), c.dot(a), c.dot(b)])
True

If the dimension of your vectors is higher, then you can pick any vector that is not in the plane defined by a,b and subtract that projection, the normalize.

Example in higher dimensions

# 1. setup
n = 5
np.random.seed(0)
a, b = np.random.uniform(size=(2, n))
a /= np.linalg.norm(a)
b -= a.dot(b)*a
b /= np.linalg.norm(b)

# 2. check
assert np.allclose([1,1,0,0], [a.dot(a), b.dot(b), a.dot(b), b.dot(a)])

Then:

# 3. solve
ab = np.c_[a, b]
c = np.roll(a   b, 1)  # any vector unlikely to be 0 or some
                       # linear combination of a and b
c -= (c @ ab) @ ab.T
c /= np.linalg.norm(c)

# 4. check
abc = np.c_[a, b, c]
>>> np.allclose(np.eye(3), abc.T @ abc)
True

Generalization: complement an m-basis in a n-D space

In an n-dimensional space, given an (n, m) orthonormal basis x with m s.t. 1 <= m < n (in other words, m vectors in a n-dimensional space put together as columns of x): find n - m vectors that are orthonormal, and that are all orthogonal to x.

We can do this in one shot using SVD.

# 1. setup
#    we use SVD for the setup as well, for convenience,
#    but it's not necessary at all. It is sufficient that
#    x.T @ x == I

n, m = 6, 2  # for example
x, _, _ = np.linalg.svd(np.random.uniform(size=(n, m)))
x = x[:, :m]

# 2. check
>>> np.allclose(x.T @ x, np.eye(m))
True

>>> x.shape
(6, 2)

So, at this point, x is orthonormal and of shape (n, m).

Find y to be one (of possibly many) orthonormal basis that is orthogonal to x:

# 3. solve
u, s, v = np.linalg.svd(x)
y = u[:, m:]

# 4. check
>>> np.allclose(y.T @ y, np.eye(n-m))
True

>>> np.allclose(x.T @ y, 0)
True