Multilinear operators

[mateuszbaran]

I've started some work to solve https://github.com/JuliaManifolds/ManifoldsBase.jl/issues/45 . It's far from complete but the general idea is to have a AbstractTensorField type that represent arbitrary partial functions between vector bundles, which covers gradients, Hessians, metric fields, etc.

TODO:

• [x] Metric field apply_operator.
• [ ] Metric field get_coordinates.
• [x] Cotangent vector field apply_operator.
• [ ] Cotangent vector field get_coordinates.
• [ ] Naming things in a way that feels right.
• [ ] More extensive testing for MetricManifold.
• [ ] Connections as an additional example.
• [x] ~~Do we consider complex vector bundles on real manifolds and similar mixed scalar field cases? Probably VectorSpaceType should hold the scalar field over which the vector space is over.~~

Stretch goals:

• [ ] Raising/lowering indices by moving vector spaces between VSIn and VSOut. Includes extending sharp and flat.
• [ ] Implementing get_coordinates for operators. Takes bases of both in and out space. How would tangent space bases relate to this? I'm not sure yet.
• [ ] Jacobian get_coordinates.
• [ ] Hessian get_coordinates.
• [ ] Lowering/raising indices for Hessian example.
• [ ] Consider whether get_vector analogue makes sense for these (even if it does I would call it out-of-scope for this issue).
• [ ] Jacobian apply_operator.
• [ ] Hessian apply_operator.

[mateuszbaran]

@kellertuer could you take a look if my approach to scalars for vector spaces makes sense? There may be some holes related to scalars of tangent spaces of complex manifolds. I'm also not sure if the relation to complex bases is clear enough (it's described in one of the docstrings).

[kellertuer]

Can you point me to the code lines? I don't yet see through the new structure you are building here.

[mateuszbaran]

Some code lines:

• Here is the new definition of VectorSpaceType: https://github.com/JuliaManifolds/Manifolds.jl/pull/202/files#diff-c6389631d82343df20214336402283efR26 . It is now parametrized by number system. I replaced the constants TangentSpace and CotangentSpace by functions since they depend on the number system. Although I'm still not sure how to relate number system of a basis with a number system of a tangent space.
• I was trying to clear some things about number systems here: https://github.com/JuliaManifolds/Manifolds.jl/pull/202/files#diff-c6389631d82343df20214336402283efR62-R80 , where abstract tensor fields are described. Also here, in the docstring for TangentBundleFibers: https://github.com/JuliaManifolds/Manifolds.jl/pull/202/files#diff-1eb1c4a47b912d3f0697cda39383efcdR21-R38 .

[kellertuer]

Looks fine to me and seems to be as generic as possible, which I like :)

[mateuszbaran]

Now after some discussion in https://discourse.julialang.org/t/ann-tensoralgebra-jl-taking-covariance-seriously/43587 I'm starting to think that perhaps adding a number system argument to vector spaces wasn't a good idea. After all it can be deduced from a basis.

[kellertuer]

You're right, it can be deduced from a basis.

[mateuszbaran]

I've read some parts of the Absil's book and reevaluated (again) the plans here. It looks like the most important missing component are affine connections and that's what I've started working on in the last commit. I'm thinking about opening a separate PR with a narrower scope, with only affine connections and things that can be easily built on top of them, just like I did with gradients and differentials of curves.

[codecov[bot]]

Codecov Report

Merging #202 into master will decrease coverage by 0.58%. The diff coverage is 55.26%.

@@            Coverage Diff             @@
##           master     #202      +/-   ##
==========================================
- Coverage   95.03%   94.45%   -0.59%     
==========================================
  Files          63       62       -1     
  Lines        3744     3712      -32     
==========================================
- Hits         3558     3506      -52     
- Misses        186      206      +20     

Impacted Files	Coverage Δ
src/Manifolds.jl	100.00% <ø> (ø)
src/differentiation.jl	89.47% <0.00%> (-10.53%)	:arrow_down:
src/finite_diff.jl	83.33% <0.00%> (-16.67%)	:arrow_down:
src/vector_space.jl	33.33% <33.33%> (ø)
src/manifolds/Circle.jl	99.23% <100.00%> (ø)
src/manifolds/Euclidean.jl	92.13% <100.00%> (+0.37%)	:arrow_up:
src/manifolds/PowerManifold.jl	95.83% <100.00%> (-0.05%)	:arrow_down:
src/manifolds/ProductManifold.jl	89.12% <100.00%> (-0.45%)	:arrow_down:
src/manifolds/VectorBundle.jl	97.40% <100.00%> (-0.03%)	:arrow_down:
src/tests/tests_general.jl	98.54% <100.00%> (ø)
... and 8 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 785b1f8...4533539. Read the comment docs.

[kellertuer]

Could you also shortly check why Julia 1 is failing? It worked on my last commit before I merged to master but on master it failed; I think they might have introduced a rounding error in the complex case maybe?

[mateuszbaran]

I've increased exp_log_atol_multiplier for complex Grassmann manifold, maybe it will help.

[kellertuer]

It seems that the GitHub action is broken, at least in my log file, the Julia v1 action loads Julia 1.5.0.

[mateuszbaran]

That's weird. The CI file has 1.0 explicitly. I've added 1.5 and we'll see what happens.

And it looks like I've been fixing the Grassmann problem wrong then.

[kellertuer]

I think your approach was good, I also only noticed by chance, that the Julia 1.0 test was actually running 1.5 (and can't see why it does)

[mateuszbaran]

I'll ask about it in the github-actions channel on Slack.

And my approach was incorrect because I thought the failure is on 1.0, not 1.5... it's also on Windows so I've replaced the version check with system check.

[mateuszbaran]

github-actions channel was very fast. Version 1.0 was parsed as 1 which is then converted to the latest stable version of Julia 1 😕 . Apparently we didn't have CI for Julia 1.0 for quite some time.

[kellertuer]

Oh. I wasn't aware of that. Interesting. Well, let's hope 1.0 does not produce too many errors now.

[kellertuer]

Maybe we should fix the 1.0 issue in a separate PR and for simplicity you disable it here, too?

see #209.

[mateuszbaran]

Thanks, good idea.

[kellertuer]

I think this might help a lot, getting the Quasi newton more performant, as far as I see my student struggled there a little bit, will discuss with him this week. There is a branch but not yet a PR for (L-)BFGS in Manopt.

[mateuszbaran]

What exactly is the problem this would help with? I stopped working on it because I didn't really have any applications that could guide the design. The other thing is that it's a very complex problem. I think it would be reasonable to first do some specific operators like cotangent vectors or maybe connections and then generalize it.

By the way, do you currently have any form of ehess2rhess in Manopt.jl? Maybe that's what is needed?

[kellertuer]

While we have a way to get from egrad to rgrad using project, I think we could have a use case to get from ehess to an rhess but currently I do not have good idea how to do that (also haven'T thought about it much).

The problem we might have is, that we currently have two ways representing a hessian

• a set of tangent vectors (might require the large dimension of the embedding)
• in coordinates (requires parallel transpiring a basis) and we currently do the second approach (I can surely provide the links to code lines) and we might use too much memory. To provide context: on the sphere we beat the Matlab variant. On Stiefel(1000,2) we don't (by far, say a factor of 300).

That's why I am interested in an efficient way to do the Hessian :)

[mateuszbaran]

While we have a way to get from egrad to rgrad using project, I think we could have a use case to get from ehess to an rhess but currently I do not have good idea how to do that (also haven'T thought about it much).

Weingarten maps look like a reasonable solution. It's described here: https://link.springer.com/chapter/10.1007/978-3-642-40020-9_39 .

The problem we might have is, that we currently have two ways representing a hessian

• a set of tangent vectors (might require the large dimension of the embedding)
• in coordinates (requires parallel transpiring a basis) and we currently do the second approach (I can surely provide the links to code lines) and we might use too much memory. To provide context: on the sphere we beat the Matlab variant. On Stiefel(1000,2) we don't (by far, say a factor of 300).

Well, I think the problem is that you don't really want to represent the Hessian itself, just calculate its values on tangent vectors. On (abstract) embedded manifolds you can just calculate Euclidean Hessian-vector product using an explicit formula or an AD tool like ReverseDiff.jl and then "Riemannize" it using a Weingarten map like in that link. I haven't thought about all the details here but I hope the general idea is clear enough.

[kellertuer]

Yes, sure I want to store (as manopt does, too) the action on tangents, i.e. H maps a tangent to a tangent vector. Weingarten might be too complicated for that.

And yes, your idea is clear.

[mateuszbaran]

Why would Weingarten be too complicated for that?

[kellertuer]

Sure, in general it is not, but maybe in special cases there are better options.

[mateuszbaran]

I think there isn't anything particularly valuable in this PR and other approaches are generally more preferable.