check_spherical_jacob_tangent.py

# Geometric intepretation of the gradient of the mapping:
# f : (0, inf) x Sphere(k-1) -> R^k
# f(r, u) -> r*u
# The *catch*: R can vary on (0, inf) but u may only vary on the
# k-1--dimensional tangent plane!

import numpy as np

def main():

    k = 4
    r = .42

    rng = np.random.RandomState(42)
    x = rng.randn(k)
    x /= np.linalg.norm(x)

    # identify a basis of the tangent space to the sphere at x.
    # This gives us the valid directions of variation.
    u, s, v  = np.linalg.svd(np.eye(k) - np.outer(x, x))
    B = u.T[:-1]

    # we have basis @ x = 0

    # gradient in euclidean coordinates:
    d_eucl = r * np.eye(4)

    # gradient expressed wrt tangent basis:
    # (must solve B.T @ d_tg = d_eucl; multiply by B on both sides.)
    d_tg = B @ d_eucl
    print(d_tg.shape)
    print(d_tg.round(2))
    print(np.linalg.svd(d_tg)[1])

    # print(np.dot(d_tg, d_tg.T).round(2))
    # = diag(r ** 2).

    # the riemannian gradient should be d_tg multiplied by
    # the metric tensor at x. But **i think** because we are using
    # the ambient R^k coordinates, the metric tensor is Id.

    J = np.row_stack([x, d_tg])
    u, s, v = np.linalg.svd(J)
    print(s)
    print(u.round(2))
    print(np.linalg.det(J))
    print("expected", r ** 3)

    # can we prove this? Probably first "deflate" x which has singular value 1,
    # then show that the remaining stuff (exactly d_tg) has k-1 times the singular value r.


if __name__ == '__main__':
    main()