frechet_mean_projective_space.py

num_points = 10000
num_dims = 200

def norm(x, keepdims=False):
    return np.linalg.norm(x, ord=2, axis=-1, keepdims=keepdims)

def normalize(x):
    return x / norm(x, True)

points = normalize(np.exp(np.random.randn(num_points, num_dims)))


def loss(x, points):
    x = normalize(x)
    return np.sum(np.inner(x, points)**2, axis=-1)

def minimum(points):
    A = np.matmul(points.T, points)
    w, v = np.linalg.eigh(A)
    ind1 = np.all(v > 0, axis=0)
    ind2 = np.all(v < 0, axis=0)
    if np.any(ind1):
        ind = ind1
    elif np.any(ind2):
        ind = ind2
    else:
        raise RuntimeError
    return w[ind], np.abs(v[:, ind])

w, v = minimum(points)
print(w)#, v)
print(loss(v.T, points))
print(loss(points, points).max())