Stiefel manifold restricted to a subspace: Hessian calculation

[rs1909]

I am trying to work out how to implement a version of the Stiefel manifold, where the matrix is restricted to a subspace. In terms of an implicit equation, this would be defined by

p^T p - I = 0
B^T p = 0,

where B is a constant orthogonal matrix. I can make the projection work so P_p X = (I-B B^T)X - (p^T X + X^T p)/2
and it looks like the Polar retraction is actually the same as for the standard Stiefel.

The gradient does work, i.e., G = P_p [D f(p)], but the Hessian is a problem. Could anyone help with this? This is more of a theoretical question, because once I have the formula I can implement it.
Thanks...

[kellertuer]

HI,
I (and we in the package) do not have worked much with a Hessian yet. So for example in the differentiation section of this package there is no support for Hessians yet – just because in general they require the derivative of a vector field on the tangent bundle.

Can you elaborate what you mean with “the gradient works”? And where does the Hessian then “not work”?

[rs1909]

I am implementing a Newton's method to find an optimum of a function with such a restricted Stiefel manifold. So I have an objective function 'f', that takes an input from the set of matrices, which should actually be on this restricted Stiefel manifold. So what I mean is that I can accurately calculate the Riemannian gradient of 'f', because that is just the projection of the function derivative 'Df' to the tangent space T_p M.

I was able to work out the Hessian calculation with the standard Stiefel manifold, but for some reason when I do a similar calculation it just does not work with the restricted Stiefel, as in the Hessian is not accurate.

In Manifolds.jl you are also calculating vector transports and those have something to do with Hessians. However I never understood how to used them to get the Riemannian Hessian from D^2 f using vector transports. I think this is what Manopt.jl uses...

In any case I need to study this a bit more. I was just a bit desperate with this.

[kellertuer]

You could use a quasi-Newton approach, e.g. L-RBFGS, see for example https://manoptjl.org/stable/solvers/quasi_Newton.html – there the Hessian is approximated from the gradient.

[mateuszbaran]

You can take a look at this paper: https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.867.9051&rep=rep1&type=pdf , it describes an approach to calculation of Riemannian Hessians.

[rs1909]

Thanks, I will look into it.