CODEWITHSUNDEEP
pythonnumpylinear-algebra
a_guest
a_guest

Reputation: 36339

Generate random vector which is perpendicular to a set of given vectors

I have a set of vectors which are pairwise perpendicular. Then I would like to generate a new vector which should be perpendicular to all the vectors in the given set and it should be chosen at random. Basically, if d is the dimensionality of the problem and m is the number of given vectors, then the set of vectors that are perpendicular to those spans a (d-m)-dimensional subspace. I would like to randomly sample a vector from this subspace.

I can use np.linalg.lstsq in order to find a vector that is perpendicular to all the others, but the result won't be random:

import numpy as np

rng = np.random.default_rng(seed=0)

d = 8  # number of dimensions

vectors = [rng.normal(size=d)]
for _ in range(d-1):
    vectors.append(np.linalg.lstsq(
        np.stack(vectors, axis=0),
        np.zeros(len(vectors)),
        rcond=None,
    )[0])


import itertools as it

for i, j in it.combinations(range(d), r=2):
    assert abs(vectors[i] @ vectors[j]) < 1e-16

So I think I should be sampling (d-m) components at random and then determine the others with np.linalg.solve. For example:

import numpy as np

rng = np.random.default_rng(seed=0)

d = 10  # number of dimensions

vectors = [rng.normal(size=d)]
for _ in range(d-1):
    random_indices = rng.choice(np.arange(d), size=d-len(vectors), replace=False)
    random_mask = np.zeros(d, dtype=bool)
    random_mask[random_indices] = True

    random_components = rng.normal(size=random_mask.sum())

    matrix = np.stack(vectors, axis=0)
    other_components = np.linalg.solve(
        matrix[:, ~random_mask],
        -(matrix[:, random_mask] @ random_components),
    )
    new = np.empty_like(vectors[-1])
    new[random_mask] = random_components
    new[~random_mask] = other_components

    for v in vectors:
        assert abs(v @ new) < 1e-12, f'dot product: {v @ new}'
    vectors.append(new)

For the above to work, however, I have to relax the perpendicularity condition v @ new == 0 depending on the number of dimensions. For example for d = 10 I can only require v @ new < 1e-12 while for d = 20 the threshold is 1e-10. For the least squares solution, it was sufficient to use the same threshold independent of d (in fact, all dot products were zero).

Also, I'm not absolutely certain if the above algorithm, i.e. first randomly sampling the indices of to-be-randomized components, then determining the others, will yield a truly random vector out of this subspace. How could I verify this?

Upvotes: 1

Views: 890

Answers (1)

Bob
Bob

Reputation: 14664

You can do this iteratively using the Gram-Schmidt process

n = 100;
d = 20;
v = np.random.rand(n);
v = v / np.sqrt(np.sum(np.abs(v)**2));
V = [v]

for i in range(d):
    v = np.random.rand(n);
    # orthogonalize
    v = v - sum(vi * np.sum(vi.conj() * v) for vi in V);
    # normalize
    v = v / np.sqrt(np.sum(np.abs(v)**2))
    V.append(v);
A = np.array(V)
np.allclose(A @ A.T, np.eye(d+1)) # check the result

In this example I started with an empty basis.

If you already have the basis and want a single vector you simply repeat what I did in the last iteration.

Upvotes: 0

Related Questions