`interior_point_Newton` on `ProductManifold`

[JoshuaLampert]

I'm trying to solve a constrained optimization problem with the interior_point_Newton method on a ProductManifold and come across two problems (so far), where I think I was able to solve the first one, but am still struggling with the second one.
The MWE I am using is:

using Manifolds, Manopt
using RecursiveArrayTools

N = 5
M1 = SkewSymmetricMatrices(N)
M2 = PositiveVectors(N)
M = ProductManifold(M1, M2)

# Just some dummy objective and constraint
f(M, x) = 1.0
grad_f(M, x) = 0.0
h(M, x) = 0.0
grad_h(M, x) = 0.0
objective = ConstrainedManifoldObjective(f, grad_f; h = h, grad_h = grad_h, g = nothing,
                                         grad_g = nothing, equality_constraints = 1)

p = rand(M)
interior_point_Newton(M, objective, p)

For the first one I needed to adjust Manopt.jl in three places, namely:

Manopt.jl/src/plans/interior_point_Newton_plan.jl

Line 1015 in 44e2bfc

N = M × vector_space(length(μ)) × vector_space(length(λ)) × vector_space(length(s)),

Manopt.jl/src/solvers/interior_point_Newton.jl

Line 190 in 44e2bfc

_step_M::AbstractManifold = M × vector_space(length(μ)) × vector_space(length(λ)) ×

and

Manopt.jl/src/solvers/interior_point_Newton.jl

Line 207 in 44e2bfc

_sub_M = M × vector_space(length(λ)),

where in each place I needed to replace the × by explicit calls to ProductManifold (e.g., for the first one N = ProductManifold(M, vector_space(length(μ)), vector_space(length(λ)), vector_space(length(s))). As I understand this is necessary for M being already a ProductManifold because the code relies on things like N[1] being M again, which in the × case would only return the first component of M and not the ProductManifold itself. However, by using ProductManifold(M, ...) one gets a nested ProductManifold, which seems to me as the wanted behavior in this case. If you agree, I am happy to create a PR for that.
After resolving this, I get the stacktrace

ERROR: LoadError: MethodError: no method matching get_hessian(::ProductManifold{…}, ::ManifoldFirstOrderObjective{…}, ::ArrayPartition{…}, ::ArrayPartition{…})
The function `get_hessian` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  get_hessian(::AbstractManifold, ::EmbeddedManifoldObjective{P, Missing, E} where E, ::Any, ::Any) where P
   @ Manopt ~/.julia/dev/Manopt/src/plans/embedded_objective.jl:139
  get_hessian(::AbstractManifold, ::EmbeddedManifoldObjective{P, T, E} where E, ::Any, ::Any) where {P, T}
   @ Manopt ~/.julia/dev/Manopt/src/plans/embedded_objective.jl:151
  get_hessian(::AbstractManifold, ::ScaledManifoldObjective, ::Any, ::Any)
   @ Manopt ~/.julia/dev/Manopt/src/plans/scaled_objective.jl:105
  ...

Stacktrace:
 [1] get_hessian(M::ProductManifold{…}, co::ConstrainedManifoldObjective{…}, p::ArrayPartition{…}, X::ArrayPartition{…})
   @ Manopt ~/.julia/dev/Manopt/src/plans/constrained_plan.jl:805
 [2] (::CondensedKKTVectorFieldJacobian{…})(N::ProductManifold{…}, Y::ArrayPartition{…}, q::ArrayPartition{…}, X::ArrayPartition{…})
   @ Manopt ~/.julia/dev/Manopt/src/plans/interior_point_Newton_plan.jl:432
 [3] (::CondensedKKTVectorFieldJacobian{…})(N::ProductManifold{…}, q::ArrayPartition{…}, X::ArrayPartition{…})
   @ Manopt ~/.julia/dev/Manopt/src/plans/interior_point_Newton_plan.jl:420
 [4] get_hessian(TpM::Fiber{…}, slso::SymmetricLinearSystemObjective{…}, X::ArrayPartition{…}, V::ArrayPartition{…})
   @ Manopt ~/.julia/dev/Manopt/src/plans/conjugate_residual_plan.jl:149
 [5] get_gradient(TpM::Fiber{…}, slso::SymmetricLinearSystemObjective{…}, X::ArrayPartition{…})
   @ Manopt ~/.julia/dev/Manopt/src/plans/conjugate_residual_plan.jl:113
 [6] interior_point_Newton(M::ProductManifold{…}, cmo::ConstrainedManifoldObjective{…}, p::ArrayPartition{…}; kwargs::@Kwargs{…})
   @ Manopt ~/.julia/dev/Manopt/src/solvers/interior_point_Newton.jl:142
 [7] top-level scope
   @ ~/.julia/dev/SummationByPartsOperatorsExtra/run/debug_Manopt.jl:89
 [8] include(mapexpr::Function, mod::Module, _path::String)
   @ Base ./Base.jl:307
 [9] top-level scope
   @ REPL[1]:1
in expression starting at /home/lampert/.julia/dev/SummationByPartsOperatorsExtra/run/debug_Manopt.jl:89
Some type information was truncated. Use `show(err)` to see complete types.

which comes from trying to create the sub_state here:

Manopt.jl/src/solvers/interior_point_Newton.jl

Lines 222 to 231 in 44e2bfc

	sub_state::St = decorate_state!(
	ConjugateResidualState(
	TangentSpace(_sub_M, _sub_p),
	sub_objective;
	X = _sub_X,
	stop = sub_stopping_criterion,
	sub_kwargs...,
	);
	sub_kwargs...,
	),

Reducing this even further one can also get the same error by calling

Xp = rand(M)
get_hessian(M, objective, p, Xp)

instead of the whole call to interior_point_Newton in the MWE above. Is there a missing definition of get_hessian or do you have an idea how to fix this?

[kellertuer]

For the first part: Yes that sounds reasonable, in short because × “appends” to a product manifold, so if one of the factors is a product manifold, the new one is a product manifold of one element longer. Since that is not what we want here, yes, we have to write that as an explicit ProductManifold which does not do the unwrapping.
That both behave differently is on purpose, since both might be a behaviour someone needs; but at time of writing we did not have that in mind the correct way around,

[kellertuer]

The second part I do not yet understand, and while we had exciting stuff we worked on this week, most of that was coding and debugging, so my brain is a bit fried to hunt bugs for the next few days.

At first glance: A firstorderobjective only has first order information, that is - by type, this does not provide second order information like a Hessian – but that for now only describes why the error appears in the last step, not yet where it comes from.

[kellertuer]

Ah!

You are entering IPNewton with only first order information (using the middle layer f directly providing the objective), not that the IPNewton needs access to a Hessian, see the docs at

https://manoptjl.org/stable/solvers/interior_point_Newton/

which states

or a ConstrainedManifoldObjective cmo containing f, grad_f, Hess_f, and the constraints

It is a bit hard to constrain the type in the objective layer that you use (generating the objective before and passing that directly), that is why the high-level function does not error (one could maybe add checks), but you need a Hessian or an approximate one at least, see https://manoptjl.org/stable/solvers/trust_regions/#Manopt.ApproxHessianBFGS for one approximate hessian that we do offer.

So – ys for the first please do a PR.
For the second we should maybe look how to improve the error message and error earlier, if no Hessian access is available; you can for now resolve it by either implementing a Hessian or using the approx above.

edit: You also need Hess_h in your case. Since which hessians you need depends on the constraints, that is also a bit hard to check on the “we pass an objective” level. But if you have ideas how to improve the error message, that is also welcome :)

[kellertuer]

Ah and yes the math is broken in https://manoptjl.org/stable/solvers/interior_point_Newton/#Manopt.interior_point_Newton as well, I will try to find time to fix that unless you feel adventurous to fix the math as well ;)

[JoshuaLampert]

Thanks for your quick reply! I'll try to pass a Hessian and see if that works. I also opened #541 for the first point.

[kellertuer]

Nice, thanks! Will take a look tomorrow.

Keep in mind that the algorithm and its implementation here are both relatively new (both 08/24); I still hope it is stable enough already.

[JoshuaLampert]

Alright. Good to know. And no worries. I'll see if it satisfies my needs and maybe make the implementation a bit more stable along the way :)

[kellertuer]

Your PR is just on the new version to be registered right now, so you can use it directly soon.

Let me know when you encounter more things or whether this issue can then be closed.

[JoshuaLampert]

Nice, thanks! Yes, I will try explicitly passing the Hessian later today and will report back.

[JoshuaLampert]

We are getting closer. I now run into the following four different problems/questions.

1. I modified my initial dummy problem to

using Manifolds, Manopt, ForwardDiff, ADTypes
using RecursiveArrayTools

N = 5
M1 = SkewSymmetricMatrices(N)
M2 = PositiveVectors(N)
M = ProductManifold(M1, M2)

function optimization_gradient(M, f, p, autodiff)
    b = Manifolds.ManifoldDiff.TangentDiffBackend(autodiff)
    return Manifolds.ManifoldDiff.gradient(M, p -> f(M, p), p, b)
end

autodiff = AutoForwardDiff()
f(M, p) = 1.0
grad_f(M, p) = optimization_gradient(M, f, p, autodiff)
hess_f(M, p, Xp) = ApproxHessianBFGS(M, p, grad_f)(M, p, Xp)
h(M, p) = 0.0
grad_h(M, p) = optimization_gradient(M, h, p, autodiff)
hess_h(M, p, Xp) = ApproxHessianBFGS(M, p, grad_h)(M, p, Xp)
objective = ConstrainedManifoldObjective(f, grad_f; hess_f = hess_f,
                                         h = h, grad_h = grad_h, hess_h = hess_h,
                                         g = nothing, grad_g = nothing,
                                         equality_constraints = 1)

p = rand(M)
x = interior_point_Newton(M, objective, p)

where I already use AD for the gradient because this is what I will need later for my real problem. This now at least runs through, but I get only NaNs as solution. Using the debug keyword, shows that the :GradientNorm is NaN already starting in the first iteration. I would expect it to be zero. Am I missing something or is this a bug?
2. This is probably a general question about equality constraints and not specifically related to using ProductManifolds with interior_point_Newton. In my true problem I get the following error

ERROR: LoadError: The point ArrayPartition{Float64, Tuple{Matrix{Float64}, Vector{Float64}}}(([0.0 0.6756502488724243 -0.2666666666666667 0.14101641779424268 -0.05; -0.6756502488724243 0.0 0.950460144138989 -0.41582631306080786 0.14101641779424268; 0.2666666666666667 -0.950460144138989 0.0 0.9504601441389892 -0.2666666666666667; -0.14101641779424268 0.41582631306080786 -0.9504601441389892 0.0 0.6756502488724242; 0.05 -0.14101641779424268 0.2666666666666667 -0.6756502488724242 0.0], [0.1, 0.5444444444444446, 0.7111111111111111, 0.5444444444444446, 0.1])) on ProductManifold(SkewSymmetricMatrices(5), PositiveVectors(5)) is not feasible for the provided constants.

* There are 0 of 0 inequality constraints violated. 
* There are 1 of 1 equality constraints violated. The sum of violation is 4.437342591868191e-31.

On the one hand this makes sense because 4.437342591868191e-31 is not zero and therefore the constraint is numerically violated. On the other hand this is a bit frustrating because 4.437342591868191e-31 is obviously very close to 0.0 and since we are comparing floating point numbers we usually don't expect exact equality anyway. Since my constraint is a bit more difficult, it would be extremely hard for me to find an initial guess, which exactly satisfies the constraint without rounding error. Would it make sense to introduce a tolerance when checking for feasibility or suppress any feasibility checking for the initial guess? After commenting out

Manopt.jl/src/solvers/interior_point_Newton.jl

Line 243 in c10dec8

!is_feasible(M, cmo, p; error = :error)

and

Manopt.jl/src/solvers/interior_point_Newton.jl

Line 265 in c10dec8

!is_feasible(M, cmo, ips.p; error = :error)

the code runs (but I got the NaNs again (starting in the second iteration); maybe due to the same reason as in point 1?).
3. If I let the optimization for my real problem run for only one iteration to inspect the non-NaN result after the first iteration, I see that the equality constraint is not satisfied as h(M, p) gives 0.23832525161749957, where p is the solution. Is it expected that results can violate the constraint or would you say we cannot trust the solution anyway since it will turn NaN in the next iteration? I understand it might be difficult to tell without the code, but I am not (yet) ready to share my real code and am more interested in your general assessment about solutions violating equality constraints.
4. This is not so important when the approximate Hessian already works, but I was wondering if also using AD for the Hessians work. I tried to use ManifoldDiff.hessian with the TangentDiffBackend, but in the dummy example above something like

b = Manifolds.ManifoldDiff.TangentDiffBackend(autodiff)
Manifolds.ManifoldDiff.hessian(M, p -> f(M, p), rand(M), b)

gives me a (zero) 15 by 15 matrix. I am not sure how to interpret this. I would have expected something like an ArrayPartition of two N = 5 by N = 5 matrices, which can be multiplied component-wise to an element Xp of M, which itself is an ArrayPartition of an N by N matrix and a vector of length N such that the result is again an ArrayPartition of an N by N matrix and a vector of length N .

[kellertuer]

1. I do not yet have any idea, since I am not using AD tools much (nearly none at all), so that would probably be more of a deep dive for me once I find time.

2. I think I do not fully understand where the error is from, since you did not provide a stacktrace, but any violation below 1.0e-16 I would consider to be a rounding error.

3. no an interior point method stays interior (feasible) but also only if you start in the interior. So you really have to make sure your start point is feasible. More than that I can not say in this general setup.

4. While I have used an AD gradient a bit, I have never used a Hessian, maybe Check with Mateusz by opening an Issue on ManifoldDiff. In you example I am not sure what b is for example? If that is a basis because you want a Hessian matrix (and not like usually in Manopt the Hessian in the sense of the Hessian with a direction – or in the Euclidean case H*v instead of the Hessian matrix H).
   For a Hessian matrix you always need a basis for the tangent space (consisting of as many basis vectors as the dimension of the manifold) and the Hessian matrix will always be with respect to this basis (and coordinate vectors w.r.t. that, because otherwise this is not representable.

As an example, for a function f on the sphere S2 we do represent points and tangent vectors as 3D vectors, but all problems we work on are of cause 2D.
Any tangent vector can hence be represented instead of by X (in R3) also by coefficients c /in R2) as soon as we fix a basis in a tangent space (say b1, b2) by decomposing X = c1b1 + c2b2. This is a one-to-one correspondence.
Now grad f(p) makes sense in both X and c, but it is easier to work with in X
For the Hessian we usually implement Y = (Hess f(p))[X], which you can think of as X -> c then a matrix multiplication of the Hessian in coordinates d = Hc, and then d -> Y, just that the direct way is usually much easier to do.
But it is not easy to trite the direct way (for app p and X) as an actual easy matrix multiplication. For the actual metric you need the tangent space basis. In that basis only the whole length makes sense.

I do not yet see your relation from 15 to 5 and 5, but I probably miss a few details on your example there.

[kellertuer]

More details on this, best to read Section 5 of Nicolas book https://www.nicolasboumal.net/#book as well as 8.11 for the generic / abstract case.

[JoshuaLampert]

First of all, thanks a lot for the detailed and quick answer!
2. The error comes from the two lines I quoted (that's why commenting them out fixes the issue) calling

Manopt.jl/src/plans/constrained_plan.jl

Line 1004 in c10dec8

function is_feasible(M, cmo, p; check_point::Bool = true, error::Symbol = :none, kwargs...)

Just to be clear: I do not want to understand why this error is thrown. I am very well aware of why (which is why I didn't post a complete stacktrace). It's because feasibility of equality constraints are compared for exact equality here (probably should have made that more clear before, sorry):

Manopt.jl/src/plans/constrained_plan.jl

Line 1008 in c10dec8

feasible = v && all(g .<= 0) && all(h .== 0)

and here:

Manopt.jl/src/plans/constrained_plan.jl

Line 1041 in c10dec8

h_violated = sum(h .!= 0)

I agree with you that any violation below 1e-16 should be treated as rounding error. Therefore, my question would mostly be: Do you consider it reasonable to replace the two lines 1008 and 1041 in constrained_plan.jl I cited above by calls to isapprox with tolerances instead of testing for equality. If yes, I am happy to provide a PR for that.
3. Yes, that's what I would have expected, but it's not what is happening here. My initial point is feasible (well, if a value for h of 4.437342591868191e-31 counts as feasible, which we agree on, see point 2), but the result is not. That's why I was confused and asked. But probably something funky is going on anyway because of the NaNs quickly appearing. So before thinking more about this point, one would need to fix that first.
4. I'll stick with the approximate Hessian for now. That's good enough for me right now. Again, the NaN issue should be solved first. To answer your question

In you example I am not sure what b is for example?

I am not sure if that is what you mean, but my b is defined in the code snippet:

b = Manifolds.ManifoldDiff.TangentDiffBackend(autodiff)

It's just the backend, which is needed by ManifoldDiff.hessian.

[kellertuer]

For 2,
one could pass down / keep / model some tolerances (the question is then where to store them, but sure that is something we can check. Also the is_point would accept several kwargs (like an atol).

We can check the next days how to best do that; not sure when exactly I find time.

For 3, I would need more input, out of the blue I have no clue.

For 4, ah then I was not reading careful enough, I was expecting some basis to be passed as well. Maybe easiest to do an MWE and post a separate issue to ManifoldDiff for that one, still.