Robust Procrustes using standard (manifold unaware) Levenberg-Marquardt vs "Riemannian" Levenberg-Marquardt

[ojwoodford]

@mateuszbaran @kellertuer You recently published a paper on "A MODIFIED RIEMANNIAN LEVENBERG-MARQUARDT ALGORITHM FOR ROBUST OR CONSTRAINT OPTIMIZATION ON MANIFOLDS"

Some optimizers, such as L-BFGS, transport gradient information from one point in the state space to another, so they need to be modified to work on manifolds. However, other optimizers, including Levenberg-Marquardt, do not; they only use gradient information computed at the current state, and use it to update the current state. It's therefore easy to hide the existence of the manifold from the optimizer, using a "box-plus" interface and thus use an off the shelf optimizer. My question is, why wouldn't you do this? How does the optimizer (in this case Levenberg-Marquardt) benefit from knowing about the manifold?

To make the point, I implemented the Robust Procrustes problem from here, in an open-source, manifold agnostic, implementation of Levenberg-Marquardt from the NLLSsolver package. Here are graphs of the result (with code here).

The results show that the two implementations converge to the same result, and that the standard implementation converges in far fewer iterations. I believe these results are equivalent to what Ceres (which you cited but didn't compare against) would have achieved. The runtime isn't so useful, because it mainly depends on the speed of the Jacobian computation, which can of course be improved if necessary. My code uses AD, which transitions from statically to dynamically sized arrays, at the point where my method becomes slower.

I suppose this is about two things. The first is not reinventing the wheel. If standard Levenberg-Marquardt works fine, why claim to have modified it, so it now works? Why not benefit from existing software that already works well. The second is about APIs. Manifolds is a useful package. Why not make it compatible with a "box plus" interface?

I had a not very fruitful discussion with @kellertuer about the box-plus interface. I understood that he doesn't understand exactly what it is, but also thinks it is not very good, and he doesn't want to discuss it further. This struck me as an odd position for an academic working in this area to take, given that it works perfectly well (better even than your solution). I'm happy to talk more about it, but only if people are interested. Otherwise I'll leave it there.

[kellertuer]

It is an interesting aspect, that you claim that I do not understand box-plus. I said I heavily dislike it. And I mentioned all arguments already, mainly its limited use.

Storing gradients is not enough – “moving on the manifold” is the other important part. Plot+project (what box plus does) does not always work, is probably very inexact.
Manifolds.jl has therefore usually several retractions – and when available the exponential map, to check which one works best.
And I wrote that before. That is the reason I dislike box plus. But ti seems I have to write it again.

For the examples I wrote – Sphere, Robust procrustes, Robust Subspace procrustes.
The Sphere is merely an introuctionary example that aims to compare to the non robust case. There is no need to compare there.
Neither SO(d), d > 3 not Stiefel are available in Ceres. Also for SO(3) they use a completely different representation.
So I am very unsure how Ceres should have obtained them.
At first glance you use Static arrays. We indeed did not. That usually has an
Since your code is merely undocumented, I am not so sure what it does.

We do not claim to reinvent LM, we generalise the robust case to manifolds.

You seem to be very convinced your code is super much better. What do you want to achieve? That I stop developing Manopt?
Yes, I considered that a few times, since feedback like yours I get regularly – and you bring me again closer to that. Is that what you intend to aim for here? After all these repetitions of a nearly annoying level, I do not see any other goal you aim to achieve.

So please stop doing that and also stop to claim things you do just assume.

[kellertuer]

I'm very happy to talk more about it with anyone who's interested.

From the very beginning your tone arrives to me as "I can do everything better, why even bother doing yours? You should use mine!"
And that is too aggressive for me to be interested in a discussion, especially when I am still convinced, yours is a very special case of ours - we could do that (given more time) - as well as to me, still, yours intention is totally unclear.

Take a look, check whether Manopt achieves the same with static arrays and projection, if not, sure then you are indeed faster in a special case.

Still unclear to me what the goal then is (besides aggressively proposing we stop our approach).

I do not want to continue in such an approach that from the beginning had this aggressive tone, so this is indeed my last comment here.

[ojwoodford]

I would say I found the method to work worse than standard LM - more iterations, and more inconsistency in the final cost. You didn't have a proper baseline in your paper either for that example. I'm pointing these things out. It's constructive criticism. Take it or leave it.

[ojwoodford]

Take a look, check whether Manopt achieves the same with static arrays and projection, if not, sure then you are indeed faster in a special case.

I didn't make any point about the speed. I deliberately said it wasn't relevant.

[ojwoodford]

I do not want to continue in such an approach that from the beginning had this aggressive tone, so this is indeed my last comment here.

I think my tone has generally been far less aggresive than yours, Ronny.

[ojwoodford]

aggressively proposing we stop our approach

I haven't said that. As I said before, I encourage you to stop making judgements about what I think or am trying to achieve. If you want to engage, engage with what I actually say.

[mateuszbaran]

It's therefore easy to hide the existence of the manifold from the optimizer, using a "box-plus" interface and thus use an off the shelf optimizer. My question is, why wouldn't you do this? How does the optimizer (in this case Levenberg-Marquardt) benefit from knowing about the manifold?

Why would we hide the manifold from the solver if it can benefit from it by for example constraining the step size by injectivity radius? What do we get by hiding part of the problem from the solver? Sure, the implementation is slightly more simple but not that much, and some solvers have Euclidean-only tricks but these are quite rare among nonlinear solvers.

How does the optimizer (in this case Levenberg-Marquardt) benefit from knowing about the manifold?

For LM it's mostly stepsize limiting and being able to support arbitrarily complex representation of iterates and their updates in a way that is transparent to the solver.

. It's therefore easy to hide the existence of the manifold from the optimizer, using a "box-plus" interface and thus use an off the shelf optimizer

Box-plus is effectively just a retraction. It can't work with off-the-shelf optimizers that are not manifold-aware. Using box-plus interface is a basic example of manifold-awareness.

The results show that the two implementations converge to the same result, and that the standard implementation converges in far fewer iterations. I believe these results are equivalent to what Ceres (which you cited but didn't compare against) would have achieved.

Robust Procrustes is a simple problem, and LM is a fairly simple algorithm. If two different implementations take different numbers of iterations, it's most likely due to differences in hyperparameters.

We don't have a comparison with Ceres because we haven't managed to construct a meaningful examples where both libraries can be easily used and we have a clear advantage: either we roughly match their performance, they have an advantage due to being hyper-focused on a use case, or we have a feature that they don't like the matrix-free solving of linear problems or bounds support. I personally didn't spend time looking for such example because this solver is far enough from my primary research that I didn't need it for anything particular.

The first is not reinventing the wheel. If standard Levenberg-Marquardt works fine, why claim to have modified it, so it now works? Why not benefit from existing software that already works well. The second is about APIs. Manifolds is a useful package. Why not make it compatible with a "box plus" interface?

Please let me know which library supports all the features our LM does, including matrix-free solves and bounds. I haven't found one, and if there isn't one -- no wheel was reinvented. There are problems for which both features become relevant.

The box-plus interface is a subset of Manifolds.jl interface. There is no incompatibility.

I had a not very fruitful discussion with @kellertuer about the box-plus interface. I understood that he doesn't understand exactly what it is, but also thinks it is not very good, and he doesn't want to discuss it further. This struck me as an odd position for an academic working in this area to take, given that it works perfectly well (better even than your solution). I'm happy to talk more about it, but only if people are interested. Otherwise I'll leave it there.

We know what that is, it's just that we are already familiar with a much more rich interface and don't see a good reason to restrict ourselves to this subset chosen specifically for one solver.

[ojwoodford]

Box-plus is effectively just a retraction. It can't work with off-the-shelf optimizers that are not manifold-aware. Using box-plus interface is a basic example of manifold-awareness.

Box-plus is not defined to be anything particular - it can be a retraction or exponentiation. A box-plus interface just overloads the plus operator. It's the kind of composition that Julia excels at, and the optimizer is absolutely unaware of the manifold. This is how my code works.

[mateuszbaran]

Any box-plus operation can be represented by our retract.

absolutely unaware of the manifold

Supplying a solver with box-plus makes it indirectly aware of the manifold. This is a pointless distinction that makes it more difficult to coordinate multiple manifold-related operations.

[ojwoodford]

Supplying a solver with box-plus makes it indirectly aware of the manifold. This is a pointless distinction that makes it more difficult to coordinate multiple manifold-related operations.

We'll have to disagree on the optimizer's awareness of the manifold. I think my point is about composability. A standard optimizer can be used with variables on a manifold if the interface to the state space allows the overloading of the plus operator. The whole value in Julia is the composability. More composability means fewer implementations, easier maintenance, and more contributors.

Please let me know which library supports all the features our LM does, including matrix-free solves and bounds.

The paper was about Riemannian LM with robustness and bounds. Ceres does that, notwithstanding the implementation of the n-d rotation manifold, which is straightforward.

You say that in fact the optimizer can benefit from knowledge of the manifold. I'm genuinely interested to know how. You've stated the following:

matrix-free solves

I don't know about this. Can you give an example? Was such an example used in the paper?

makes it more difficult to coordinate multiple manifold-related operations

I'm not sure what you mean. Can you give an example?

constraining the step size by injectivity radius

Is this necessary? Does it improve things? The step size is already constrained to reduce the cost.

If two different implementations take different numbers of iterations, it's most likely due to differences in hyperparameters

It's certainly possible, but "most likely", I disagree on. It's interesting to note that there is a clear pattern to the number of iterations and the final cost to my solver, while yours is far more variable. That suggests a more fundamental issue to me. For comparison, below are graphs using the dogleg optimizer rather than LM, again using default hyperparameters.

[ojwoodford]

Also from Claude: "Can ceres-solver support bounds and manifolds on the same problem?"

Yes. In Ceres, bounds and manifolds are both per-parameter-block properties, and nothing stops you setting both on the same problem — or even the same block. Internally, ParameterBlock::Plus applies Manifold::Plus(x, delta) first, then clamps the result onto the box constraints (the code comment says exactly: "Manifold::Plus() followed by projection onto the hyper cube implied by the bounds constraints").
Two caveats:

1. The common, safe case is bounds on some blocks and manifolds on others — that works fine (with a trust-region minimizer; the line-search minimizer doesn't support bounds).
2. Combining both on the same block is mathematically dubious for something like Stiefel: the clamp happens in ambient coordinates after the retraction, so it will generally knock the point off the manifold (clamping entries of an orthonormal matrix breaks orthonormality), and the clamping isn't reflected in the Jacobian. It only makes sense when the box projection maps the manifold to itself. If you need genuine constraints on top of a Stiefel block, the usual tricks are reparameterization or penalty terms instead.

[mateuszbaran]

Yes — the best match is MatrixManifoldsInCeres by Tolga Birdal.

Thanks, I didn't now about that library. Still, it's just a subset of what we can do.

[ojwoodford]

I'm not questioning the usefulness of Manifolds.jl at all. I'm questioning whether every optimizer ever written (or even just the good ones) needs to be "adapted" to work on manifolds, or whether, as I currently believe, they should work out of the box. From my brief testing, it seems the canonical (standard) optimizers work better.

[mateuszbaran]

Yes. In Ceres, bounds and manifolds are both per-parameter-block properties, and nothing stops you setting both on the same problem — or even the same block. Internally, ParameterBlock::Plus applies Manifold::Plus(x, delta) first, then clamps the result onto the box constraints (the code comment says exactly: "Manifold::Plus() followed by projection onto the hyper cube implied by the bounds constraints").

I didn't see it anywhere in the docs. Anyway, we do GCP/L-BFGS-B-style search on the model instead of a simple projection.

[mateuszbaran]

I'm not questioning the usefulness of Manifolds.jl at all. I'm questioning whether every optimizer ever written (or even just the good ones) needs to be "adapted" to work on manifolds, or whether, as I currently believe, they should work out of the box. From my brief testing, it seems the canonical (standard) optimizers work better.

Ceres will be better for >95% of problems than our implementation. So most people don't need it. I'm personally not quite sure the work on RLM was worth it, though it's done and might be useful to someone at some point. I'm not planning on working on adding manifolds to any more solvers if that makes you feel better.

[ojwoodford]

It occurred to me that a difference in number of iterations might be due to the way that the update is parameterized. In my code the update is parameterized according to this line of code:

NLLSsolver.update(var::Rotation{N, N2, T}, updatevec, start=1) where {N, N2, T} = Rotation(var.R * proj2rot(vec2skew(updatevec, start, N)))

The current NxN rotation matrix is var.R, and updatevec is an N(N-1)/2 long update vector. The update vector is turned into a first order approximation of a small rotation by vec2skew, by creating a skew-symmetric matrix and adding the identity. That matrix is then projected into the space of rotation matrices by proj2rot. Finally the current rotation is composed with update rotation.

How is the update parameterized in your approach?

This highlights a benefit of a box-plus approach: It's easy see, and to change, the way an update is parameterized, and this can be done on a per variable block basis. For example, I could choose to retract the state rotation, rather than just the update rotation, as follows:

NLLSsolver.update(var::Rotation{N, N2, T}, updatevec, start=1) where {N, N2, T} = Rotation(proj2rot(var.R * vec2skew(updatevec, start, N)))

I expect that would perform worse.

[mateuszbaran]

This is all possible in our approach as well, and much more. Please read the paper and docs before guessing and making questions that are already answered.

[ojwoodford]

Actually, it didn't perform any differently.

[ojwoodford]

This is all possible in our approach as well, and much more. Please read the paper and docs before guessing and making questions that are already answered.

I can see that LevenbergMarquardt takes a retraction_method keyword argument, but it isn't clear to me how to define different retraction methods for different variables within the optimization.