Adding RDoG

[Nimrais]

I recently played around in my internal lab with the RDoG optimizer https://arxiv.org/pdf/2406.02296 (more methods are proposed in the paper I just do not have time to try them all), which is basically a Riemannian generalization of the DoG optimizer https://arxiv.org/abs/2302.12022. First, it's very easy to understand, and second, it's not computationally intensive compared to Armijo. My cost functions are usually very expensive to evaluate—surprisingly, sometimes gradients are comparatively cheap (at least the stochastic gradients)—so it would be nice to have something that is not prohibitively slow but at the same time learning-rate schedule free. If you are finding it valuable to have it in the repo I can file a PR.

My current implementation is the following one:

@doc raw"""
    RDoGStepsize{R<:Real} <: Stepsize
"""
mutable struct RDoGStepsize{R<:Real,P} <: Stepsize
    initial_distance::R
    max_distance::R
    gradient_sum::R
    initial_point::P
    use_curvature::Bool
    sectional_curvature_bound::R
    last_stepsize::R
end

function RDoGStepsize(
    M::AbstractManifold;
    p=rand(M),
    initial_distance::R=1e-3,
    use_curvature::Bool=false,
    sectional_curvature_bound::R=0.0,
) where {R<:Real}
    return RDoGStepsize{R,typeof(p)}(
        initial_distance,
        initial_distance,  # max_distance starts at initial_distance
        zero(R),          # gradient_sum starts at 0
        copy(M, p),       # store initial point
        use_curvature,
        sectional_curvature_bound,
        NaN,              # last_stepsize
    )
end

"""
    geometric_curvature_function(κ::Real, d::Real)
"""
function geometric_curvature_function(κ::Real, d::Real)
    if κ < 0 && d > 0
        sqrt_abs_κ = sqrt(abs(κ))
        arg = sqrt_abs_κ * d
        # Avoid numerical issues for small arguments
        if arg < 1e-8
            return 1.0 + arg^2 / 3.0  # Taylor expansion
        else
            return arg / tanh(arg)
        end
    else
        return 1.0
    end
end

function (rdog::RDoGStepsize)(
    mp::AbstractManoptProblem,
    s::AbstractManoptSolverState,
    i,
    args...;
    gradient=nothing,
    kwargs...,
)
    M = get_manifold(mp)
    p = get_iterate(s)
    grad = isnothing(gradient) ? get_gradient(mp, p) : gradient
    
    # Compute gradient norm
    grad_norm_sq = norm(M, p, grad)^2
    
    if i == 0
        # Initialize on first call
        rdog.gradient_sum = grad_norm_sq
        rdog.initial_point = copy(M, p)
        rdog.max_distance = rdog.initial_distance
        
        # Initial stepsize
        if rdog.use_curvature
            ζ = geometric_curvature_function(rdog.sectional_curvature_bound, rdog.max_distance)
            stepsize = rdog.initial_distance / (sqrt(ζ) * sqrt(max(grad_norm_sq, eps())))
        else
            stepsize = rdog.initial_distance / sqrt(max(grad_norm_sq, eps()))
        end
    else
        # Update gradient sum
        rdog.gradient_sum += grad_norm_sq
        
        # Update max distance
        current_distance = distance(M, rdog.initial_point, p)
        rdog.max_distance = max(rdog.max_distance, current_distance)
        
        # Compute stepsize
        if rdog.use_curvature
            ζ = geometric_curvature_function(rdog.sectional_curvature_bound, rdog.max_distance)
            stepsize = rdog.max_distance / (sqrt(ζ) * sqrt(rdog.gradient_sum))
        else
            stepsize = rdog.max_distance / sqrt(rdog.gradient_sum)
        end
    end
    
    rdog.last_stepsize = stepsize
    return stepsize
end

get_initial_stepsize(rdog::RDoGStepsize) = rdog.last_stepsize
get_last_stepsize(rdog::RDoGStepsize) = rdog.last_stepsize

function show(io::IO, rdog::RDoGStepsize)
    s = """
    RDoGStepsize(;
      initial_distance=$(rdog.initial_distance),
      use_curvature=$(rdog.use_curvature),
      sectional_curvature_bound=$(rdog.sectional_curvature_bound)
    )
    
    Current state:
      max_distance = $(rdog.max_distance)
      gradient_sum = $(rdog.gradient_sum)
      last_stepsize = $(rdog.last_stepsize)
    """
    return print(io, s)
end

"""
    RDoG(; kwargs...)
    RDoG(M::AbstractManifold; kwargs...)

Creates a factory for the [`RDoGStepsize`](@ref).

## Keyword arguments

* `initial_distance=1e-3`: initial distance estimate
* `use_curvature=false`: whether to use sectional curvature
* `sectional_curvature_bound=0.0`: lower bound on sectional curvature

$(_note(:ManifoldDefaultFactory, "RDoGStepsize"))
"""
function RDoG(args...; kwargs...)
    return ManifoldDefaultsFactory(Manopt.RDoGStepsize, args...; kwargs...)
end

I forked the repo and made an example for the Rayleigh quotient https://github.com/Nimrais/Manopt.jl/blob/master/examples/rdog_example.jl.

If you will call it you will obtain the following result.

True optimal value (negative largest eigenvalue): -103.56226794615476
RDoG Example: Rayleigh Quotient Minimization
==================================================
Manifold: Sphere(29)
True optimal value: -103.56226794615476
Initial cost: -26.050530058561083
Initial optimality gap: 77.51173788759368

1. Basic RDoG (no curvature)
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -108.764589219706 | ||grad f|| = 5.46144063e+00 | step = 7.14544209e-03
# 100    | f(x) = -111.794090201414 | ||grad f|| = 1.94207092e+01 | step = 4.83908048e-03
# 150    | f(x) = -105.322098782001 | ||grad f|| = 3.73581804e+00 | step = 4.75430550e-03
# 200    | f(x) = -103.914644314998 | ||grad f|| = 7.31942441e-01 | step = 4.75104968e-03
# 250    | f(x) = -103.631999324778 | ||grad f|| = 1.44217419e-01 | step = 4.75092292e-03
# 300    | f(x) = -103.576035552189 | ||grad f|| = 2.84495965e-02 | step = 4.75091798e-03
# 350    | f(x) = -103.564984971476 | ||grad f|| = 5.61355494e-03 | step = 4.75091779e-03
# 400    | f(x) = -103.562804101333 | ||grad f|| = 1.10769530e-03 | step = 4.75091778e-03
# 450    | f(x) = -103.562373744711 | ||grad f|| = 2.18578122e-04 | step = 4.75091778e-03
# 500    | f(x) = -103.562288823125 | ||grad f|| = 4.31314309e-05 | step = 4.75091778e-03
# 550    | f(x) = -103.562272065753 | ||grad f|| = 8.51101171e-06 | step = 4.75091778e-03
# 600    | f(x) = -103.562268759064 | ||grad f|| = 1.67945553e-06 | step = 4.75091778e-03
Final cost: -103.56226842977263
Final optimality gap: -4.836178675304836e-7
Gap reduction factor: 1.6027476049095297e8

2. RDoG with curvature
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -72.095369960602 | ||grad f|| = 3.70251339e+01 | step = 5.96437367e-03
# 100    | f(x) = -99.259701666648 | ||grad f|| = 8.38849476e+00 | step = 4.75172205e-03
# 150    | f(x) = -103.000821847999 | ||grad f|| = 1.15867074e+00 | step = 4.73360379e-03
# 200    | f(x) = -103.486796766971 | ||grad f|| = 1.56925069e-01 | step = 4.73327489e-03
# 250    | f(x) = -103.552080484853 | ||grad f|| = 2.12039156e-02 | step = 4.73326889e-03
# 300    | f(x) = -103.560892014525 | ||grad f|| = 2.86421892e-03 | step = 4.73326878e-03
# 350    | f(x) = -103.562082096827 | ||grad f|| = 3.86881848e-04 | step = 4.73326878e-03
# 400    | f(x) = -103.562242842924 | ||grad f|| = 5.22574324e-05 | step = 4.73326878e-03
# 450    | f(x) = -103.562264555380 | ||grad f|| = 7.05858209e-06 | step = 4.73326878e-03
Final cost: -103.5622684228676
Final optimality gap: -4.767128416460764e-7
Gap reduction factor: 1.6259628672881514e8

3. Fixed stepsize for comparison
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -12.164350964363 | ||grad f|| = 6.14247524e+01 | step = 1.00000000e-02
# 100    | f(x) = -154.102934108739 | ||grad f|| = 6.36359713e+01 | step = 1.00000000e-02
# 150    | f(x) = -27.826871769201 | ||grad f|| = 1.21849235e+02 | step = 1.00000000e-02
# 200    | f(x) = -142.833567182374 | ||grad f|| = 7.66632630e+01 | step = 1.00000000e-02
# 250    | f(x) = -11.489914079089 | ||grad f|| = 6.21234688e+01 | step = 1.00000000e-02
# 300    | f(x) = -73.826859643045 | ||grad f|| = 6.10347964e+01 | step = 1.00000000e-02
# 350    | f(x) = -47.379687572027 | ||grad f|| = 4.76704604e+01 | step = 1.00000000e-02
# 400    | f(x) = -129.266585427988 | ||grad f|| = 1.80502748e+01 | step = 1.00000000e-02
# 450    | f(x) = -145.390758234062 | ||grad f|| = 7.51024945e+01 | step = 1.00000000e-02
# 500    | f(x) = -32.820973565456 | ||grad f|| = 1.24955303e+02 | step = 1.00000000e-02
# 550    | f(x) = -90.420049394265 | ||grad f|| = 8.14205164e+00 | step = 1.00000000e-02
# 600    | f(x) = -143.388053266504 | ||grad f|| = 7.63607133e+01 | step = 1.00000000e-02
# 650    | f(x) = -77.428020925510 | ||grad f|| = 1.59325019e+01 | step = 1.00000000e-02
# 700    | f(x) = -13.069167334947 | ||grad f|| = 6.04488467e+01 | step = 1.00000000e-02
# 750    | f(x) = -15.525608524218 | ||grad f|| = 1.12794817e+02 | step = 1.00000000e-02
# 800    | f(x) = -138.759517588186 | ||grad f|| = 2.63777228e+01 | step = 1.00000000e-02
# 850    | f(x) = -71.881686168002 | ||grad f|| = 6.01177083e+01 | step = 1.00000000e-02
# 900    | f(x) = -133.778383487263 | ||grad f|| = 2.18157716e+01 | step = 1.00000000e-02
# 950    | f(x) = -125.908532239189 | ||grad f|| = 7.95445936e+01 | step = 1.00000000e-02
# 1000   | f(x) = -114.380405799096 | ||grad f|| = 7.72801332e+01 | step = 1.00000000e-02
Final cost: -114.38040579909566
Final optimality gap: -10.818137852940893
Gap reduction factor: 7.164979679614846

4. Sensitivity to initial distance
------------------------------
Running shorter tests (200 iterations) with different initial distances:
ε = 1.0e-5: Final gap = 0.2792590185977133, Reduction = 277.56216532169816x
ε = 0.001: Final gap = -0.037547838633443575, Reduction = 2064.3461969753585x
ε = 0.1: Final gap = 0.30871421932334897, Reduction = 251.07926048073432x

[kellertuer]

I think every step size idea is valuable and can be added :)

For your code I only have a few minor comments:

• sure we should extend the doc strings a bit and include some formulae
• One could maybe use a bit longer name that conveys a bit more information? Maybe like DistanceOverGradients?
• For the show method, I split that into a show for just the constructor and a status_summary which is printed e.g. in the show of a solver status and would also report the internal values you state.

But all the rest looks great! You even already included the Factory!

Independent of the question to add it (if you find time we surely should add it!)
I am not so much a fan of a comparison to a constant step size, unless you optimise that to really fine tune the constant and find nearly the best constant step size. How does it compare to Armijo?

So yes, feel free to do a PR, remember to also do

• tests
• add yourself to the about.md that lists all contributions (that are more than fixing typos)
• feel free to add yourself to the zenodo.json as well then.

edit: You are right the rule is not that complicated – even the weighted one looks reasonably easy to do. So we could even do both. For the constant step size – for example their example in the paper is also just horrible. When people come along and only know about constant step size and are angry about how bad it performs – well maybe do some math ;) ? There is no guarantees for that, for Armijo there are! But maybe that is also just my usual rambling about ML people when they forget that math does exist ;)

[Nimrais]

Hi! Thanks for the quick and thoughtful feedback — happy to open a PR.

I’ll expand the docstrings and include the formulae. I can also rename the type to something more descriptive, e.g. DistanceOverGradientsStepsize, while keeping RDoG as the factory for convenience (or providing both names)

I will split the solver printing into show(::RDoGStepsize) (constructor-style) and status_summary(::RDoGStepsize) so the solver can display internal state cleanly.

You’re absolutely right that comparing against a single constant stepsize is fragile unless optimized per problem. I’ve added an Armijo comparison on the sphere in the example so it’s apples-to-apples. And Armijo generaly does better job. But I must to say that Armijo as an algorithm is in another computational budget category. Armijo uses a backtracking routine that repeatedly evaluates the cost during the increase/decrease loops, for my experiments it's just not vary practical. In contrast, RDoG’s stepsize uses only the gradient norm and distance from the initial point — no extra cost evaluations — which is why it’s attractive when gradients are cheaper than costs.

This gives a clearer picture: Armijo tends to reduce the gap more aggressively but at higher per-iteration cost due to backtracking evaluations, whereas RDoG is learning-rate free and very cheap per iteration. However, in this comparison below, RDoG beats Armijo. Another issue is that it's sometimes not very practical to ask about sectional curvature as well. However, I'm always questioning my hyperparameter choices for Armijo, so it could be just a skill issue on my side and someone could do far better. But that's also kind of a fruitful aspect of RDoG—that you need only one hyperparameter.

True optimal value (negative largest eigenvalue): -103.56226794615476
RDoG Example: Rayleigh Quotient Minimization
==================================================
Manifold: Sphere(29)
True optimal value: -103.56226794615476
Initial cost: -26.050530058561083
Initial optimality gap: 77.51173788759368

1. Basic RDoG (no curvature)
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -108.764589219706 | ||grad f|| = 5.46144063e+00 | step = 7.14544209e-03
# 100    | f(x) = -111.794090201414 | ||grad f|| = 1.94207092e+01 | step = 4.83908048e-03
# 150    | f(x) = -105.322098782001 | ||grad f|| = 3.73581804e+00 | step = 4.75430550e-03
# 200    | f(x) = -103.914644314998 | ||grad f|| = 7.31942441e-01 | step = 4.75104968e-03
# 250    | f(x) = -103.631999324778 | ||grad f|| = 1.44217419e-01 | step = 4.75092292e-03
# 300    | f(x) = -103.576035552189 | ||grad f|| = 2.84495965e-02 | step = 4.75091798e-03
# 350    | f(x) = -103.564984971476 | ||grad f|| = 5.61355494e-03 | step = 4.75091779e-03
# 400    | f(x) = -103.562804101333 | ||grad f|| = 1.10769530e-03 | step = 4.75091778e-03
# 450    | f(x) = -103.562373744711 | ||grad f|| = 2.18578122e-04 | step = 4.75091778e-03
# 500    | f(x) = -103.562288823125 | ||grad f|| = 4.31314309e-05 | step = 4.75091778e-03
# 550    | f(x) = -103.562272065753 | ||grad f|| = 8.51101171e-06 | step = 4.75091778e-03
# 600    | f(x) = -103.562268759064 | ||grad f|| = 1.67945553e-06 | step = 4.75091778e-03
Final cost: -103.56226842977263
Final optimality gap: -4.836178675304836e-7
Gap reduction factor: 1.6027476049095297e8

2. RDoG with curvature
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -72.095369960602 | ||grad f|| = 3.70251339e+01 | step = 5.96437367e-03
# 100    | f(x) = -99.259701666648 | ||grad f|| = 8.38849476e+00 | step = 4.75172205e-03
# 150    | f(x) = -103.000821847999 | ||grad f|| = 1.15867074e+00 | step = 4.73360379e-03
# 200    | f(x) = -103.486796766971 | ||grad f|| = 1.56925069e-01 | step = 4.73327489e-03
# 250    | f(x) = -103.552080484853 | ||grad f|| = 2.12039156e-02 | step = 4.73326889e-03
# 300    | f(x) = -103.560892014525 | ||grad f|| = 2.86421892e-03 | step = 4.73326878e-03
# 350    | f(x) = -103.562082096827 | ||grad f|| = 3.86881848e-04 | step = 4.73326878e-03
# 400    | f(x) = -103.562242842924 | ||grad f|| = 5.22574324e-05 | step = 4.73326878e-03
# 450    | f(x) = -103.562264555380 | ||grad f|| = 7.05858209e-06 | step = 4.73326878e-03
Final cost: -103.5622684228676
Final optimality gap: -4.767128416460764e-7
Gap reduction factor: 1.6259628672881514e8

3. Fixed stepsize for comparison
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -12.164350964363 | ||grad f|| = 6.14247524e+01 | step = 1.00000000e-02
# 100    | f(x) = -154.102934108739 | ||grad f|| = 6.36359713e+01 | step = 1.00000000e-02
# 150    | f(x) = -27.826871769201 | ||grad f|| = 1.21849235e+02 | step = 1.00000000e-02
# 200    | f(x) = -142.833567182374 | ||grad f|| = 7.66632630e+01 | step = 1.00000000e-02
# 250    | f(x) = -11.489914079089 | ||grad f|| = 6.21234688e+01 | step = 1.00000000e-02
# 300    | f(x) = -73.826859643045 | ||grad f|| = 6.10347964e+01 | step = 1.00000000e-02
# 350    | f(x) = -47.379687572027 | ||grad f|| = 4.76704604e+01 | step = 1.00000000e-02
# 400    | f(x) = -129.266585427988 | ||grad f|| = 1.80502748e+01 | step = 1.00000000e-02
# 450    | f(x) = -145.390758234062 | ||grad f|| = 7.51024945e+01 | step = 1.00000000e-02
# 500    | f(x) = -32.820973565456 | ||grad f|| = 1.24955303e+02 | step = 1.00000000e-02
# 550    | f(x) = -90.420049394265 | ||grad f|| = 8.14205164e+00 | step = 1.00000000e-02
# 600    | f(x) = -143.388053266504 | ||grad f|| = 7.63607133e+01 | step = 1.00000000e-02
# 650    | f(x) = -77.428020925510 | ||grad f|| = 1.59325019e+01 | step = 1.00000000e-02
# 700    | f(x) = -13.069167334947 | ||grad f|| = 6.04488467e+01 | step = 1.00000000e-02
# 750    | f(x) = -15.525608524218 | ||grad f|| = 1.12794817e+02 | step = 1.00000000e-02
# 800    | f(x) = -138.759517588186 | ||grad f|| = 2.63777228e+01 | step = 1.00000000e-02
# 850    | f(x) = -71.881686168002 | ||grad f|| = 6.01177083e+01 | step = 1.00000000e-02
# 900    | f(x) = -133.778383487263 | ||grad f|| = 2.18157716e+01 | step = 1.00000000e-02
# 950    | f(x) = -125.908532239189 | ||grad f|| = 7.95445936e+01 | step = 1.00000000e-02
# 1000   | f(x) = -114.380405799096 | ||grad f|| = 7.72801332e+01 | step = 1.00000000e-02
Final cost: -114.38040579909566
Final optimality gap: -10.818137852940893
Gap reduction factor: 7.164979679614846

4. Sensitivity to initial distance
------------------------------
Running shorter tests (200 iterations) with different initial distances:
ε = 1.0e-5: Final gap = 0.2792590185977133, Reduction = 277.56216532169816x
ε = 0.001: Final gap = -0.037547838633443575, Reduction = 2064.3461969753585x
ε = 0.1: Final gap = 0.30871421932334897, Reduction = 251.07926048073432x

5. Armijo (line search) for comparison
------------------------------
Initial  | f(x) = -26.050530058561 |  | # 50     | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 3.54960363e-17
# 100    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 150    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 200    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 250    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 300    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 350    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 400    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 450    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 500    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 550    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 600    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 650    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 700    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 750    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 800    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 850    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 900    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 950    | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
# 1000   | f(x) = -103.562236873746 | ||grad f|| = 5.68487925e-01 | step = 2.00771358e-18
Final cost (Armijo): -103.56223687374609
Final optimality gap (Armijo): 3.1072408674504004e-5
Gap reduction factor (Armijo): 2.4945519576406307e6

[kellertuer]

Armijo should have stopped loooong before that. I am surprised it does not get the gradient small – is it the right gradient (i.e. the Euclidean one)?

For consistency I would prefer the factory and the stepsize having the same name up to a Stepsize in the end for the actual. stepwise. all other step size rules follow that naming scheme as well.

[Nimrais]

Armijo should have stopped loooong before that. I am surprised it does not get the gradient small – is it the right gradient (i.e. the Euclidean one)?

I do not see mistake in my gradient sorry, can you point me it out, maybe I am doing smt stupid?

f(x) = - p' A p then gradient -2Ap no? The gradients are also ortogonal to the p so it seems my riemannian gradient is also correct.

# Define the cost function (negative Rayleigh quotient)
f(M, p) = -p' * A * p

# Define the Riemannian gradient
function grad_f(M, p)
    euc_grad = -2 * A * p
    return euc_grad - dot(euc_grad, p) * p  # Project to tangent space
end

You can see the full example with Armijo here https://github.com/Nimrais/Manopt.jl/blob/master/examples/rdog_example.jl.

For consistency I would prefer the factory and the stepsize having the same name up to a Stepsize in the end for the actual. stepwise. all other step size rules follow that naming scheme as well.

Ah sure no problem, I will do.

[kellertuer]

hm the gradient looks right at first glance, one could check with check_gradient whether that is also stable. But well, maybe also Armijo is not that stable here then. Sometimes that happens, though I did not yet have that for Rayleigh.

[Nimrais]

But from this experiment, it seems that Armijo just has a really, really conservative estimate for the stepsize because this point does not move at all, so the gradient consequently does not decrease at all. I am not trying to generalize it beyond this specific example, however, it seems to me that on average Armijo should win.

[kellertuer]

No, the step size is numerically collapsed already before your iterate 100, which is usually an indication that your descent direction is not descend – or in other words our gradient is wrong – especially when it is still as large as in that print. On the sphere I first would expect that only for gradients less than 1e-8 in norm or so.

[Nimrais]

Okay I see, I think you are right I found some time to play with check_gradient.

using Manopt
using Manifolds
using LinearAlgebra
using Random
using ForwardDiff

# Problem: Rayleigh quotient on Sphere(d-1)
Random.seed!(42)
d = 2
M = Sphere(d - 1)

A = randn(d, d); A = A' * A  # SPD

# Cost
f(M, p) = -p' * A * p

# Gradient (original)
function grad_f_original(M, p)
    euc_grad = -2 * A * p
    return euc_grad - dot(euc_grad, p) * p
end

# Gradient (projection form; should match original)
grad_f_project(M, p) = project(M, p, -2 * A * p)

# ForwardDiff Euclidean gradient (ambient), projected to tangent
# g(x) = f(M, x) with x in R^d (no projection inside the function!)
g(x) = -x' * A * x
grad_fd_euclid(x) = ForwardDiff.gradient(g, x)
grad_f_fd_project(M, p) = project(M, p, grad_fd_euclid(p))

# Point & direction
p = rand(M)                        # random point on the sphere
X = rand(M; vector_at=p)           # random tangent direction

# Sanity: all variants should be tangent
g_orig = grad_f_original(M, p)
g_proj = grad_f_project(M, p)
g_fd   = grad_f_fd_project(M, p)

println("p ∈ M?                   ", is_point(M, p))
println("X ∈ T_pM?                ", is_vector(M, p, X))
println("Grad ∈ T_pM?             ", is_vector(M, p, g_orig))
println("||orig - proj||          = ", norm(M, p, g_orig - g_proj))
println("||orig - FD-proj||       = ", norm(M, p, g_orig - g_fd))
println("||proj - FD-proj||       = ", norm(M, p, g_proj - g_fd))
println()

# A tiny directional derivative spot-check (forward difference along retraction)
Xn = X / norm(M, p, X)
t  = 1e-6
f_diff = (f(M, retract(M, p, t * Xn)) - f(M, p)) / t
println("Directional check (rough):")
println("  inner(orig, X)         = ", inner(M, p, g_orig, Xn))
println("  inner(FD-proj, X)      = ", inner(M, p, g_fd,   Xn))
println("  (f(retr(p,tX))-f(p))/t = ", f_diff)
println()

# Gradient checks (first-order expansion):
# Use a moderately tight step window for robustness
limits = (-6.0, -1.0)
Npts   = 81

println("=== check_gradient: original, default retraction ===")
ok1 = check_gradient(
    M, f, grad_f_original, p, X;
    check_vector=true,
    error=:info,
    limits=limits,
    N=Npts,
)

println("\n=== check_gradient: projected, default retraction ===")
ok2 = check_gradient(
    M, f, grad_f_project, p, X;
    check_vector=true,
    error=:info,
    limits=limits,
    N=Npts,
)

println("\n=== check_gradient: FD-projected, default retraction ===")
ok3 = check_gradient(
    M, f, grad_f_fd_project, p, X;
    check_vector=true,
    error=:info,
    limits=limits,
    N=Npts,
)

println("\n=== check_gradient: original, exponential retraction ===")
ok4 = check_gradient(
    M, f, grad_f_original, p, X;
    check_vector=true,
    error=:info,
    limits=limits,
    N=Npts,
    retraction_method=ExponentialRetraction(),
)

println("\n=== check_gradient: FD-projected, exponential retraction ===")
ok5 = check_gradient(
    M, f, grad_f_fd_project, p, X;
    check_vector=true,
    error=:info,
    limits=limits,
    N=Npts,
    retraction_method=ExponentialRetraction(),
)

println("\nSummary:")
println("  original / default        => ", ok1)
println("  projected / default       => ", ok2)
println("  FD-projected / default    => ", ok3)
println("  original / exponential    => ", ok4)
println("  FD-projected / exponential=> ", ok5)

The result is negative, the gradients are not passing the gradient check.

p ∈ M?                   true
X ∈ T_pM?                true
Grad ∈ T_pM?             true
||orig - proj||          = 0.0
||orig - FD-proj||       = 2.7755575615628914e-17
||proj - FD-proj||       = 2.7755575615628914e-17

Directional check (rough):
  inner(orig, X)         = -0.10487556969847267
  inner(FD-proj, X)      = -0.10487556969847264
  (f(retr(p,tX))-f(p))/t = -0.1048754949550812

=== check_gradient: original, default retraction ===

=== check_gradient: projected, default retraction ===

=== check_gradient: FD-projected, default retraction ===

=== check_gradient: original, exponential retraction ===

=== check_gradient: FD-projected, exponential retraction ===

Summary:
  original / default        => false
  projected / default       => false
  FD-projected / default    => false
  original / exponential    => false
  FD-projected / exponential=> false

I am not really sure what else I should do, could it be that there is a bug in the check_gradient itself?

[kellertuer]

With plots you can also plot the error (setting plot=true. I can try to find some time to look at it myself.
check_gradient is quite stable, I am sure there is no bug in there.

[kellertuer]

I am not so sure what your code does, but if I just take your p, X and grad_f_orginal, using Plots I get

julia> check_gradient(M, f, grad_f_original, p, X; plot=true)

julia> check_gradient(M, f, grad_f_original, p, X)
true

So that looks fine with me.

[Nimrais]

Hmm I guess I misused the options of the check_gradient, thanks, I will read about them later.

Currently, it's produces me the same plot as it does for you.

using Manopt
using Manifolds
using LinearAlgebra
using Random
using Plots

# Problem: Rayleigh quotient on Sphere(d-1)
Random.seed!(42)
d = 30
M = Sphere(d - 1)

A = randn(d, d); A = A' * A  # SPD

# Cost
f(M, p) = -p' * A * p

# Gradient (original)
function grad_f_original(M, p)
    euc_grad = -2 * A * p
    return euc_grad - dot(euc_grad, p) * p
end

# Gradient (projection form; should match original)
grad_f_project(M, p) = project(M, p, -2 * A * p)

# Point & direction
p = rand(M)
X = rand(M; vector_at=p)
check_gradient(M, f, grad_f_original, p, X; plot=true)

[Nimrais]

But I think we can safely conclude that, at least in this case, the performance of the RDoG is not coming from the fact that Armijo is provided with a wrong gradient. As a sanity check, I do not specify the parameters for the Armijo right now, and it does not change the results.

[kellertuer]

Yeah I do. not clearly see why Armijo “collapses” but it usually should not.

Anyways, I also do not directly see where you would add that in a PR, so independent of this discussion feel free to provide DiftanceOverGradients and if you like DistanceOverWeightedGradients as a PR :)

[kellertuer]

Finished in #505