Divergence of a vector field

[leonardogalliano]

Dear team,
I am currently working on a project where I need to compute the divergence of a vector field $f:M\to T_xM$ for a generic manifold $M$ defined using Manifolds.jl. After reviewing your package, I found it promising, but unfortunately, I was unable to achieve my goal. Please note that I am not an expert in Riemannian geometry, so my confusion could be related to a misunderstanding on my part.

To compute the divergence, I attempted to directly trace the Jacobian (though I realize this might not be the best approach). I started with a point $x$ on a sphere and defined a function $f(x,A) = Ax$ that, by construction, always stays in the tangent space ($A$ is skew-symmetric).

In Manifolds.jl, a point on the sphere is represented in embedding coordinates (i.e., as a 3D vector) rather than in 2D spherical coordinates. Since the coordinates are in Cartesian form, I expected the Jacobian of $f$ to simply be $A$ (as $f$ is linear), resulting in a $3\times3$ matrix with trace zero (which would correspond to the divergence of $f$).

However, when I compute it using ManifoldDiff, the result is a $2\times2$ matrix. I initially thought this might be because it computes the Jacobian in intrinsic coordinates, but the trace (which gives the divergence) should not depend on the coordinates. Also, the operation seems to be significantly slower compared to directly computing the Jacobian with ForwardDiff.

Could you help clarify why this is happening? Is there a way to generalise the idea to an arbitrary manifold?

Minimal working example:

using Random, LinearAlgebra, Manifolds, ManifoldDiff, ForwardDiff, DifferentiationInterface

f(x, A) = A * x

seed = 42
rng = Xoshiro(seed)
manifold = Sphere(2)
x = rand(rng, manifold)  # 3D vector
@assert is_point(manifold, x) # x is a point on the sphere
A = [0.0 -1.0 -2.0; 1.0 0.0 -3.0; 2.0 3.0 0.0]  # 3x3 matrix
@assert A' == -A  # A is skew-symmetric: A * x preserves the norm of x
@assert is_vector(manifold, x, f(x, A)) # f(x, A) is in the tangent space

J = ManifoldDiff.jacobian(
    manifold,
    TangentSpace(manifold, x),
    x -> f(x, A),
    x,
    TangentDiffBackend(AutoForwardDiff())
)

size(J)  # (2, 2)
tr(J)    # Non-zero trace

[kellertuer]

Hi,
nice that you are using our package. First of all it is a bit dangerous to mix representation of things with its actual dimension. Computing with the sphere in the embedding (representing it with 3D vectors) just makes a lot of formulae very nice. Compared to that, in spherical coordinates a lot of formulae require special cases and modulo and such.

Still the manifold is 2-dimensional (at every point locally it is isomorphic to R2, or it “looks like” a 2D plane locally).

The first problem I see with your function $f$ is that would need a very special form for the matrix A for $f$ to be a tangent vector field. Especially in A*x you would need that it “rotates” x to be ortogonal to x (to be in the tangent space) and then apply 2 2D affine transform (and nothing else).

But let's assume we have that for the following.

Before diving into the Jacobian, let's maybe first look at the differential of $f$ that would map from TxM (the way you want to “infinitesimally. change” to TxM (since we can identify Tx(TxM) trivially with TxM), So then Df(x) : TxM -> TxM – or in other words Df(x)[V] for a tangent vector V at x is a tangent vector again. And it is a linear map. (or Df(x)[aV + bW) = aDf(x)[V] + bDf(x)[W]).

Since the Sphere is a 2dim manifold, TxM is a 2D vector space. So if we choose a basis B_1, B_2 of TxM (2 linear independent tangent vectors even in the sense of linearly independent in the embedding) – we can write any V = a_1B_1+a_2B_2

Similarly we can decompose the result Df(x)[V] with respect to this basis (or even any other basis)

That way, we can compute C_{ij} as the by Df(x)[B_i] and decomposing that into the basis and look at the coefficient upfront B_j.

In these coefficients of the basis B_1,B_2 (decompose V and Df(x)[V] w.r.t. this basis) you can represent the differential as a 2x2 matrix.
That is the Jacobian.

(in R^n take the default basis of a map F: R^n -> R^n in all this description and you get the classical jacobian) – this is the Jacobian that is computed.
The same is done on arbitrary manifolds.

[leonardogalliano]

Thanks a lot for the answer!

Yes, I see that this function requires an ad-hoc form of A, it was just an example. In the code to be sure I would do something like

f_tangent(x, p) = project(manifold, x, f(x, p))

and work with f_tangent only.

I see why J is 2D then, very clear. But then, why do I get an "unexpected" divergence? I've realised that it is probably not the right way of computing it, but that I should use

$$ \mathrm{div}\left(f(x,A)\right) = \frac{1}{\sqrt{\left|G(x)\right|}}\mathrm{Tr}\left(\frac{\partial\sqrt{\left|G(x)\right|}f\left(x,A\right)}{\partial x}\right), $$

where I guess the derivatives are in intrinsic coordinates. Is there an easy way to implement this?

[kellertuer]

The remark on $A$ was mainly for completeness, good that you have that in mind :)

why do I get an "unexpected" divergence?

I can not fully answer this, because I am not sure what you expect so it is hard to answer where your unexpected-ness comes from. I am also not 100% sure where your formula for the divergence comes from, but I can tell you that neither differential nor jacobian are are metric dependent (maybe a bit if you take the Jacobian with respect to an ONB in TxM)

In your formula, where I am also not so sure where it is from (and there exist a bit too many different notations in differential geometry, so I try my best guess) – G(x) (my guess: the metric tensor matrix at x?) would as a matrix be basis-dependent but its determinant (| |?) not.
and the partial derivative is probably also the Jacobian thought of as in coordinates of an ONB.

The Jacobian above, as mentioned, is metric independent, so from that Jacobian to the divergence you have to respect the metric tensor as well, we have some support for that, see https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/metric/#Implementation-of-Metrics

[leonardogalliano]

Sorry for being a bit sloppy 🙈

I got that formula from this paper, that is what I'm trying to implement. The equation is number (3), part of Proposition 2. There they say that yes, G(x) is the matrix representation of the matrix tensor at x. And the partial derivative should be a matrix, at least dimensionally.

I will take a look at the link, thanks a lot :)

[kellertuer]

Yes sure the $\partial/\partial x$ is usually meant as the derivative with respect to parameters (e.g. in a chart $\varphi$ of $f(\varphi(p))$), but a nice connection here is that the coordinate derivatives of the chart are a basis of the tangent space (not necessarily an ONB in which things are nicer though).

[kellertuer]

...and no worries about the sloppiness. I seem to have guessed okay-ishly there :)

[mateuszbaran]

Hi! This is an interesting topic.

I see why J is 2D then, very clear. But then, why do I get an "unexpected" divergence?

I'm almost sure that's unexpected because our jacobian function assumes that f is a valid function with given domain and tangent space in some small neighborhood of x and your is not. f doesn't have to be differentiable in the embedding (this is a big point I usually make about on-manifold differentiation: naive embedding-based solutions assume that the function is differentiable in the embedding and usually fail when it's not; we can work with just differentiability on the manifold). Your f, however, is only valid at x: at other points it returns values from other tangent spaces so the result is useless.

There are a few ways to go about implementing formula (3). If you assume normal coordinates at $z(t)$ then IIRC all $G$-involving terms at $z(t)$ disappear. In a bit more general setting, with a non-moving chart, I think it was possible to express the correction term using something like our volume_density (I would have to check more carefully). In the full generality, you'd need to work with charts and currently we only have some basic functionality for them.

I attempted to implement Riemannian CNFs some time ago but I gave up before finishing, mostly due to hopelessly poor performance. Anyway, you might have better luck than me and I can give some tips.

[leonardogalliano]

Hi!
Wow, thanks! It's great that you implemented Riemannian CNFs as well. Thanks for the answer. I see now that this approach isn’t ideal.

If you assume normal coordinates at $z(t)$ then IIRC all $G$-involving terms at $z(t)$ disappear.

That’s what I suspected, and it’s why I was expecting a divergence of 0 in this case (basically coming just from the Cartesian Jacobian). But I wasn’t completely sure.

My (maybe a bit ambitious) goal is to implement a divergence function that works for (potentially) any manifold from Manifolds.jl. But I can see that handling different representations might make this tricky.

If you have any tips on how to approach this, I’d love to hear them!

[mateuszbaran]

Actually ManifoldDiffEq.jl is the primary result of my attempt to implement Riemannian CNFs. That's the part that I managed to get working reasonably well. About the CNFs themselves I realized at some point that they are going to be way too slow for the intended application no matter what so I haven't even finished implementing it. And IMO CNFs are only useful in a very narrow niche where you need both sampling and probability density estimation. If you don't need sampling, kernel density estimation works better and if you don't need probability density estimation, generative neural networks work great (we've ended up just using GANs in PyTorch).

I have the following tips:

1. Make sure that you can accept poor performance. I'd expect sampling to be about 100x slower than in a comparable neural network.
2. Maybe working at the origin of exponential charts would work fine for you. I'm not sure if it can be worked out but it would simplify things.
3. If the previous approach doesn't work, try simple charts first and Euclidean code. Quite often a single chart is enough and you get a somewhat-working prototype much easier.
4. If a single chart is not enough, I'd first try using a few of them and switch. We have an example of this technique in a tutorial for working in charts. You should have more luck with AD when working in a chart. You need to compute divergence yourself because we don't have tools for that yet (apart from maybe being able to get some inspiration from the existing Jacobian code and volume_density for exponential charts).

[leonardogalliano]

Thanks a lot. I tried a few things but I didn't have time to properly come up with something that works (mainly due to switching between representations). But I'll think more about it and keep you updated :)