Add scalar Laplacian

[juanitorduz]

Hi 👋 !
As discussed on Slack, adding the scalar Laplacian to the package could be interesting. As a reminder, the scalar Laplacian is a second-order differential operator which can be computed in local coordinates as

As suggested, the natural place to add the implementation is in the riemannian_metric module.

Now, the details:

• We need to compute the determinant of the metric. This should be easy. Note, however that we require the manifold to be orientable (for most common cases).
• In principle, we can compute the derivatives of the metric components similarly as in the christoffels method. We can take the derivatives of the determinant of the metric similarly.
• Space of admissible functions: As the Laplacian is a differential operator, it operates on smooth (or at least C2 functions if we do not want to enter the world of weak derivatives) functions. I don't know how to specify them. I guess the API should look something like

class RiemannianMetric(Connection, ABC):
    ...
    def laplacian(self, f: SmoothFunction) -> SmoothFunction:
        ...
        return Delta(f)

Where SmoothFunction is to be defined.

After looking into the details, adding differential operators can be very interesting but not straightforward, mainly because we would need to define the domain of the operators. Computing the Laplacian coefficients is relatively easy.

Any thoughts on how (and if) you would like to move forward? Do you think the method API suggested above makes sense?

[luisfpereira]

Thanks @juanitorduz, very nice way to start the discussion on this.

I kind of like the idea of the differential operators, but I'm afraid it goes a bit further from the current API. That being said, if we get more operators, maybe we can start some new module exploring this direction (instead of adding more methods to RiemannianMetric).

In terms of API, do we get a callable from laplacian that receives the arguments of f?

Additionally, I think we'll need to get such a callable for a particular point where to evaluate g (i.e. define the Laplacian at a given point):

class RiemannianMetric(Connection, ABC):
    # ...
    def laplacian(self, point, f: SmoothFunction) -> SmoothFunction:
        g = self.metric_matrix(point)
        # ...
        return Delta(f)

From a practical viewpoint, I think the code will be a bit slow, but it should be relatively easy to implement. See for example FisherRaoMetric, as I think it would use more and less the same tools.

Lastly, maybe we don't need to care much about the API for the smooth functions (we simply expect a callable), as long as we make it very clear in the docstrings the function must be smooth. This because I don't think there's a quick way of testing at runtime if the function is smooth, which means the user will always need to be responsible for passing proper functions.

Any thoughts on this @ninamiolane?