Curvature implementation #1606

Jules-Deschamps · 2022-07-15T17:51:20Z

Implementation of methods for metric and connection, so that riemannian metrics can return :

riemannian tensor
curvature
ricci tensor
scalar curvature
sectional curvature (was already implemented if curvature was)

Tests on shapes of curvature tensors symmetries of riemannian tensor on test files (geometry_test_cases and data_generation), and on negativity of sectional curvature of dirichlet distributions (on another script file, I did not know if I should add it in the test file).

…ress)

…tions

…backend tests

ninamiolane · 2022-07-16T00:48:41Z

@Jules-Deschamps Looking good! You also need to modify:https://github.com/geomstats/geomstats/blob/master/tests/data/connection_data.py
to solve this error appearing in the tests:

testing_data object doesn't have 'ricci_tensor_shape_test_data' function for pairing with 'test_ricci_tensor_shape'

You can run pytest tests/test_geomstats/test_connection.py on your machine directly (so that you do not wait for the GitHub actions to finish).

ninamiolane

LGTM when test's infrastructure problem is solved.

ninamiolane · 2022-07-16T00:51:11Z

geomstats/geometry/riemannian_metric.py

@@ -453,6 +453,25 @@ def normal_basis(self, basis, base_point=None):

        return gs.einsum("i, ikl->ikl", 1.0 / gs.sqrt(norms), basis)

+    def covariant_riemannian_tensor(self, base_point):
+        r"""Compute purely covariant version of Riemannian tensor at base_point.


Thanks, could you also add the formula for it, i.e. R_{ijkl} = g_{ll'}R_{ijk}^l' ?

Yes will do! As for the testing thing I thought I had resolved this but turns out I haven't yes, thanks!

luisfpereira

Thanks @Jules-Deschamps.

I'm just reviewing Connection.riemann_tensor now. Can you please take a look to my comments.

I'll come back to this PR soon.

luisfpereira · 2022-07-18T11:48:59Z

geomstats/geometry/connection.py

+            R_ijk^l = riemann_tensor[i,j,k,l]
+            Riemannian tensor curvature.
+        """
+        base_point = gs.to_ndarray(base_point, to_ndim=2)


Does this work when base_point is a matrix?

Can we somehow avoid its use? Because we create a new array every time we do it. christoffels should work with a base point or multiple base points, no?

Yes I was using it to know the number of points in base_point (to be able to permute the axes using gs.transpose), but if gs.move_axis (which I did not know existed) works I will get rid of these lines.

But if base_point is a matrix (representing only a point in the manifold), you will get wrong information, no?

Your branch is too far behind, that's why you cannot find moveaxis. Please merge master after pulling upstream.

It works apparently because there is a squeeze in christoffels. So for base_point of shape (1,dim), christoffels will return shape (dim, dim, dim) and jacobian_christoffels (dim, dim, dim, dim) because it is defined with a squeeze. So these lines are really just to get the number of points in base_point. I hope this is what you are talking about. This is a minimum reproducible example showing it works for several types of base_point (and namely a matrix representing only a point in the manifold).

import geomstats.backend as gs from geomstats.information_geometry.dirichlet import DirichletDistributions space = DirichletDistributions(3) point_1 = space.random_point() point_2 = gs.array([[1.0, 1.0, 1.0]]) point_3 = gs.array([[1.0, 2.0, 2.0], [1.0, 2.0, 1.0]]) print(space.metric.riemann_tensor(point_1)) print(space.metric.riemann_tensor(point_2)) print(space.metric.riemann_tensor(point_3))

Ok I will merge thanks!

So turns out I am keeping base_point = gs.to_ndarray(base_point, to_ndim=2) in my last commit also because I want to compute the jacobian christoffels for each of the points in base_point, which is not possible otherwise (or so far I haven't found anything).

What I mean here is related with #1607: riemann_tensor, the way it is coded now, will not work for manifolds whose default_point_type is matrix (in the same way christoffels doesn't).

I may be missing something in terms of concepts, but shouldn't everything also work for that case? (otherwise, as I told in #1607, we need to reconsider the class where to put this code).

base_point = gs.to_ndarray(base_point, to_ndim=2) can't work in the general case because you need to consider both point_types (matrix and vector). You need to do something like:

if self.default_point_type == 'vector': n_points = base_point.shape[0] if len(base_point.shape) > 1 else 1 else: n_points = base_point.shape[0] if len(base_points.shape > 2 else 1

or, even better:

from vectorization import get_n_points n_points = get_n_points(base_point, self.default_point_type)

geomstats/geometry/connection.py

alebrigant

Great addition 🎉 ! Just a few comments on the docstrings. Otherwise, it would be nice to have tests that check known values for the sectional curvature: for the hyperbolic space, sphere, and e.g. beta distributions.

alebrigant · 2022-07-18T19:47:05Z

geomstats/geometry/connection.py

+        return curvature
+
+    def riemann_tensor(self, base_point):
+        r"""Compute Riemannian tensor at base_point.


Maybe add in the docstring somewhere that the returned R_{ijk}^l verifies R(X_i, X_j)X_k = \sum_{l=1}^n R_{ijk}^l X_l where X_i is i-th basis vector of the tangent space.

geomstats/geometry/connection.py

alebrigant · 2022-07-18T20:01:13Z

geomstats/geometry/connection.py

+
+    def scalar_curvature(self, base_point):
+        r"""Compute scalar curvature at base_point.
+


Same: add formula of the scalar curvature in the docstring ?

alebrigant · 2022-07-18T20:52:36Z

geomstats/geometry/riemannian_metric.py

@@ -453,6 +453,25 @@ def normal_basis(self, basis, base_point=None):

        return gs.einsum("i, ikl->ikl", 1.0 / gs.sqrt(norms), basis)

+    def covariant_riemannian_tensor(self, base_point):
+        r"""Compute purely covariant version of Riemannian tensor at base_point.


This should be clarified by a formula. I suppose this is R_{ijks} = <R(X_i,X_j)X_k,X_s>=\sum_{l=1}^n R_{ijk}^lg_{ls} ? If so, it should be mentioned in the docstring, with X_i is the i-th basis vector of the tangent space.

Yes. I just checked because the tests were wrong but the convention used in the literature (or at least on wikipedia) is VERY weird (R_{ijkl} = <R(x_j, x_k)x_i, x_l>) so I think the previous successful tests were just a fluke. I will make a commit as soon as possible with your feedback!

You must beware that the sign convention taken in Wikipedia for the Riemann curvature tensor (https://en.wikipedia.org/wiki/Riemann_curvature_tensor) is not the same as the one taken in geomstats (see this docstring) and this has effects on the rest of the definitions. So all definitions should be taken from the same reference book, and not Wikipedia.

Good catch @alebrigant !! @Jules-Deschamps , you can even put a mention in all the docstrings that we do not use the same convention as in wikipedia, and therefore there are differences. Otherwise this might confuse a lot of people.

Looking at all the formulae, I'm not completely sure that they are all compatible, i.e. take the same conventions. The only way to avoid mistakes in my view, is to pick a reference, and to take all formulas from this reference. This reference should be added to all docstrings.

alebrigant · 2022-07-18T20:55:49Z

geomstats/geometry/riemannian_metric.py

+        metric = self.metric_matrix(base_point)
+        covariant_tensor = gs.einsum("...ij, ...klmj->...iklm", metric, riemann_tensor)
+        return covariant_tensor
+
    def sectional_curvature(self, tangent_vec_a, tangent_vec_b, base_point=None):
        r"""Compute the sectional curvature.


The wikipedia reference in the docstring should be removed because the convention is not the same one. Instead, one can put e.g. do Carmo, or some other book where the convention is the same.

Jules-Deschamps · 2022-07-19T00:29:01Z

I have added modifications following your feedbacks and have tested for dirichlet distributions (the sectional curvature is indeed everywhere negative), but for hyperbolic spaces it is more complicated because sectional and scalar curvatures require metric_matrix which has not yet been implemented. As for testing on the hypersphere, it is a bit complicated considering the extrinsic and intrinsic coordinates but I will try.

Also, I was wondering where I should add the tests and data generators so that I don't have to rewrite all the tests for every manifold? For example for dirichlet distributions, in dirichlet_data:

def ricci_tensor_shape_test_data(self):
        return self._ricci_tensor_shape_test_data(
            self.metric_args_list, self.space_list
        )
```, this looks like it could be replaced up above the architecture? Thanks

ninamiolane · 2022-07-19T03:23:02Z

Fantastic, thanks for addressing the comments. There are still some testing issues, which you can access by clicking on "details" in the github workflows above, or e.g. here:
https://github.com/geomstats/geomstats/runs/7400372430?check_suite_focus=true
for example:

    raise Exception(
E   Exception: testing_data object doesn't have 'covariant_riemann_tensor_bianchi_identity_test_data' function for pairing with 'test_covariant_riemann_tensor_bianchi_identity'

Yes, I believe we have to have this bit of redundant code for now:

def ricci_tensor_shape_test_data(self):
        return self._ricci_tensor_shape_test_data(
            self.metric_args_list, self.space_list
        )

but @LPereira95 was working on removing it so that you do not have to copy-paste it like this, in the future. @LPereira95 do you confirm?

alebrigant · 2022-07-19T07:37:37Z

To test sectional curvature more precisely than just the sign, you can check that for beta distributions, at any point gs.array([x, y]), it is given by

K_xy = gs.polygamma(2, x) * gs.polygamma(2, y) * gs.polygamma(2, x+y) * (
    gs.polygamma(1, x) / gs.polygamma(2, x) 
  + gs.polygamma(1, y) / gs.polygamma(2, y) 
  - gs.polygamma(1, x+y) / gs.polygamma(2, x+y)
) / detg ** 2

where

detg = gs.polygamma(1, x) * gs.polygamma(1, y) - gs.polygamma(1, x+y) * (gs.polygamma(1, x) + gs.polygamma(1, y))

See Appendix in https://www.sciencedirect.com/science/article/pii/S2405896321005851?via%3Dihub.

Jules-Deschamps · 2022-07-19T18:41:18Z

Thanks again!

I don't really understand how to address the issues pointed out by pytest because they are the same for every manifold (the data generation methods are not implemented without the '_' underscore whent they inherit from the tests so they don't know with what data to run the tests). This is not only for the bianchi identity but for all the new tests I have added in geometry_test_cases. So should I add the test methods without underscore for every manifold/geometry?

As for the testing of the sectional curvature on the Beta manifold, it works perfect! I have added the associated test.

luisfpereira · 2022-07-20T12:26:50Z

For the tests, I suggest for the moment that we add a skip to them (which we unskip as soon as I merge the tests PR). Nevertheless, notice that most of the metrics do not have metric_matrix implemented, which means most of new implementations will not work.

We still need to go over some parts of the new code (see above) before merging, as it may not be working for several cases.

alebrigant

Thanks for making the changes ! I still recommend further checks of the docstrings, to avoid all mistakes. See the detailed comments.

alebrigant · 2022-07-20T14:11:48Z

geomstats/geometry/connection.py

@@ -443,6 +443,9 @@ def curvature(self, tangent_vec_a, tangent_vec_b, tangent_vec_c, base_point):
        the curvature is defined by :math:`R(X,Y)Z = \nabla_{[X,Y]}Z
        - \nabla_X\nabla_Y Z + \nabla_Y\nabla_X Z`.

+        Convention used in the literature is:
+        R_{ijkl} = <R(x_j, x_k)x_i, x_l>


This should not appear here because we only have a connection and no Riemannian metric. However it is a good addition to RiemannianMetric.

alebrigant · 2022-07-20T14:19:51Z

geomstats/geometry/connection.py

+        -------
+        riemann_curvature : array-like, shape=[..., {dim, [n, m]}, {dim, [n, m]},
+                                                    {dim, [n, m]}, {dim, [n, m]}]
+            R_{ijk}^l = riemann_tensor[...,i,j,k,l]


Again, the convention used should be specified: what is the definition of R_{ijk}^l ?

If x_i denotes local coordinates and X_i= \partial / \partial x_i the i-th basis vector of the tangent space,

is the convention do Carmo's, where R(X_i, X_j)X_k = R_{ijk}^l X_l, which means R_{ijk}^l = dx^l(R(X_i, X_j)X_k)

or is it Wikipedia's convention, where R_{ijk}^l = dx^l(R(X_j, X_k)X_i) ?

alebrigant · 2022-07-20T14:48:52Z

geomstats/geometry/riemannian_metric.py

    def sectional_curvature(self, tangent_vec_a, tangent_vec_b, base_point=None):
        r"""Compute the sectional curvature.

        For two orthonormal tangent vectors :math:`x,y` at a base point,
        the sectional curvature is defined by :math:`<R(x, y)x, y> =
-        <R_x(y), y>`. For non-orthonormal vectors vectors, it is
+        <R_x(y), y>` (see do Carmo). For non-orthonormal vectors vectors, it is


Do you mean <R(x, y)x, y>?

This was from before me (sectional_curvature was already implemented but not curvature) but I think it is yes. I will unify the conventions if there are some incompatibilities.

alebrigant · 2022-07-20T15:03:51Z

geomstats/geometry/riemannian_metric.py

@@ -453,6 +453,25 @@ def normal_basis(self, basis, base_point=None):

        return gs.einsum("i, ikl->ikl", 1.0 / gs.sqrt(norms), basis)

+    def covariant_riemannian_tensor(self, base_point):
+        r"""Compute purely covariant version of Riemannian tensor at base_point.


Looking at all the formulae, I'm not completely sure that they are all compatible, i.e. take the same conventions. The only way to avoid mistakes in my view, is to pick a reference, and to take all formulas from this reference. This reference should be added to all docstrings.

Remove upper version limit on torch

…t data and pytorch backend code

Fix expm due to changes in scipy.linalg.expm

…ith s2

Remove greetings GitHub workflow that often fails

Fix bug in the autodiff of pullback metric matrix

ninamiolane · 2022-10-06T08:07:01Z

@Jules-Deschamps @alebrigant I'm tackling this now

review-notebook-app · 2022-10-06T08:21:27Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

Add curvature tensors - Follow-up from PR #1606

Jules-Deschamps and others added 8 commits June 29, 2022 16:35

curvature methods for connection

97fe1b3

update on connection

12da0e2

curvature tensors

fa46e5a

Add initial draft of numpy global dtype control (still a work in prog…

3338f6a

…ress)

Update numpy.linalg to control output dtype from scipy dependent func…

52ddb9a

…tions

Update _update_dtype decorator to take dtype pos as input

e3f13a8

Update decorators names, add new decorator to cast output and update …

0c1cbee

…backend tests

last updates for curvature

578ce12

Jules-Deschamps requested review from ninamiolane, luisfpereira and alebrigant July 15, 2022 17:51

better like this :) (tests checked for the wrong inequality)

0651a30

ninamiolane approved these changes Jul 16, 2022

View reviewed changes

luisfpereira requested changes Jul 18, 2022

View reviewed changes

Jules-Deschamps added 2 commits July 18, 2022 11:11

commits before merge

6ea7899

Merge branch 'master' into Jules-connection

3531927

alebrigant requested changes Jul 18, 2022

View reviewed changes

corrections for curvature + tests for dirichlet

f1190aa

Add AAC for Frechet Mean

faf3a5d

Add AAC for GPCA

0fc48a7

addition of sectional curvature testing for Beta distributions

77f7350

alebrigant requested changes Jul 20, 2022

View reviewed changes

luisfpereira mentioned this pull request Jul 21, 2022

Bug in christoffels for matrices manifolds #1607

Open

Add AAC for regression

f2bcfad

luisfpereira and others added 24 commits October 3, 2022 11:31

Merge pull request geomstats#1661 from lpereira95/torch

3f1d7eb

Remove upper version limit on torch

Add hessian to pytorch backend

7d063fe

Add inner product derivative with autodiff backend / pytorch

edfed05

Start implementing new inner product metric derivative with test, tes…

0b0decd

…t data and pytorch backend code

Derivative metric matrix works in pytorch

07e195a

Pass tests with christoffels in pytorch

ed58152

Add hardcoded tests on christoffels symbols for s2 in pullback

02f0a44

Add hessian + edit squeeze in autograd

b4eb6aa

Merge pull request geomstats#1660 from lpereira95/scipy_expm

a6346e4

Fix expm due to changes in scipy.linalg.expm

Add missing return in HyperboloidMetric.metric_matrix

512d3da

Pass derivative matrix s2 in autograd

32116f2

Pass derivative + parallel transport, slow in autograd

6c6458c

pass all tests in pullback in autograd, slow for parallel transport w…

b27df82

…ith s2

Pass in autograd and pytorch

7cf9c1c

pass derivative matrix in tensorflow

82ca331

Commenting very long test

7ac4c66

Remove leftover comments

85778ca

Address DS

8bb8b94

Address DS

6cb4731

Address comments

6ee2b1d

Address lint error

9d3534e

Remove greetings that often fails

cddce7d

Merge pull request geomstats#1666 from ninamiolane/rm-greetings

b247335

Remove greetings GitHub workflow that often fails

Merge pull request geomstats#1664 from ninamiolane/nina=bug

8efb53f

Fix bug in the autodiff of pullback metric matrix

Merge with master

9a45c8e

ninamiolane mentioned this pull request Oct 6, 2022

Add curvature tensors - Follow-up from PR #1606 #1668

Merged

8 tasks

ninamiolane closed this Oct 6, 2022

ninamiolane added a commit that referenced this pull request Oct 12, 2022

Merge pull request #1668 from ninamiolane/Jules-connection

9e5d2ab

Add curvature tensors - Follow-up from PR #1606

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Curvature implementation #1606

Curvature implementation #1606

Jules-Deschamps commented Jul 15, 2022

ninamiolane commented Jul 16, 2022

ninamiolane left a comment

ninamiolane Jul 16, 2022

Jules-Deschamps Jul 18, 2022

luisfpereira left a comment

luisfpereira Jul 18, 2022

Jules-Deschamps Jul 18, 2022

luisfpereira Jul 18, 2022

Jules-Deschamps Jul 18, 2022 •

edited

Jules-Deschamps Jul 19, 2022

luisfpereira Jul 20, 2022 •

edited

alebrigant left a comment

alebrigant Jul 18, 2022

alebrigant Jul 18, 2022

alebrigant Jul 18, 2022

Jules-Deschamps Jul 18, 2022

alebrigant Jul 19, 2022

ninamiolane Jul 19, 2022

alebrigant Jul 20, 2022

alebrigant Jul 18, 2022

Jules-Deschamps commented Jul 19, 2022

ninamiolane commented Jul 19, 2022

alebrigant commented Jul 19, 2022

Jules-Deschamps commented Jul 19, 2022

luisfpereira commented Jul 20, 2022

alebrigant left a comment

alebrigant Jul 20, 2022

alebrigant Jul 20, 2022

alebrigant Jul 20, 2022

Jules-Deschamps Jul 20, 2022

alebrigant Jul 20, 2022

ninamiolane commented Oct 6, 2022

review-notebook-app bot commented Oct 6, 2022


		def scalar_curvature(self, base_point):
		r"""Compute scalar curvature at base_point.

Curvature implementation #1606

Curvature implementation #1606

Conversation

Jules-Deschamps commented Jul 15, 2022

ninamiolane commented Jul 16, 2022

ninamiolane left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

luisfpereira left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Jules-Deschamps Jul 18, 2022 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

luisfpereira Jul 20, 2022 • edited

Choose a reason for hiding this comment

alebrigant left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Jules-Deschamps commented Jul 19, 2022

ninamiolane commented Jul 19, 2022

alebrigant commented Jul 19, 2022

Jules-Deschamps commented Jul 19, 2022

luisfpereira commented Jul 20, 2022

alebrigant left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ninamiolane commented Oct 6, 2022

review-notebook-app bot commented Oct 6, 2022

Jules-Deschamps Jul 18, 2022 •

edited

luisfpereira Jul 20, 2022 •

edited