The gradient is far from analytical in the autograd example

[Nimrais]

Hi! I wanted to try your repo.

However, I am unable to obtain the correct gradient using the autograd/finite gradient packages.

I used your example from https://manoptjl.org/dev/tutorials/AutomaticDifferentiation/.

What am I doing wrong?

using ForwardDiff
using StableRNGs
using LinearAlgebra
using Manopt, Manifolds
using FiniteDifferences, ManifoldDiff, ReverseDiff

n_dim = 200
A = randn(n_dim + 1, n_dim + 1)
A = Symmetric(A)
M = Sphere(n_dim);
p = zeros(n_dim + 1)

f1(p) = p' * A'p
gradf1(p) = 2 * (A * p - p * p' * A * p)
X1 = gradf1(p)

autograds = [
        ("forward", ManifoldDiff.FiniteDifferencesBackend()),
        ("finite", ManifoldDiff.ForwardDiffBackend()),
        ("reverse", ManifoldDiff.ReverseDiffBackend())
    ]

for (name, autograd) in autograds
    r_backend = ManifoldDiff.TangentDiffBackend(
        ManifoldDiff.FiniteDifferencesBackend()
    )
    gradf1_FD(p) = ManifoldDiff.gradient(M, f1, p, r_backend)
    p[1] = 1.0
    X2 = gradf1_FD(p)
    @show name, norm(M, p, X1 - X2)
end

producing

(name, norm(M, p, X1 - X2)) = ("forward", 28.467510473404776) (name, norm(M, p, X1 - X2)) = ("finite", 28.467510473404776) (name, norm(M, p, X1 - X2)) = ("reverse", 28.467510473404776)

[Nimrais]

My TOML

[deps] FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" ManifoldDiff = "af67fdf4-a580-4b9f-bbec-742ef357defd" Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"

[compat] FiniteDifferences = "~0.12.30" ForwardDiff = "~0.10.36" ManifoldDiff = "~0.3.6" Manifolds = "~0.8.74" Manopt = "~0.4.33" ReverseDiff = "~1.15.1" StableRNGs = "~1.0.0"

[kellertuer]

I am not 100% sure, but the first thing I did notice is that your p is chosen to be the zero vector so that is not of norm 1.

So if I set

p[1] = 1.0

after your code and run X1 and the last loop again I get

(name, norm(M, p, X1 - X2)) = ("forward", 2.2976761858423553e-12)
(name, norm(M, p, X1 - X2)) = ("finite", 2.2976761858423553e-12)
(name, norm(M, p, X1 - X2)) = ("reverse", 2.2976761858423553e-12)

So keep in mind that usually none of the functions on a manifold (here the sphere) has a defined behaviour if you plugin points not on the manifold, i.e. here we need norm 1 and other norms might have erratic behaviour.

You can always check that a point is on the manifold using

julia> is_point(M,p) # my modified point from above
true

which for your initial point yields

julia> is_point(M, zeros(n_dim+1))
false

and you can always get more information / an error message by passing

julia> is_point(M, zeros(n_dim+1), true)
ERROR: DomainError with 0.0:
The point [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] does not lie on the Sphere(200, ℝ) since its norm is not 1.

depending on your version – the next release 0.15 of ManifoldsBase.jl which is not yet available with Manopt.jl but will be the next few days, there is also is_point(M, zeros(n_dim+1); error=:warn) that instead prints that as a warning.

[Nimrais]

@kellertuer Yes, thank you. Your response makes total sense to me. The problem for me is that I am computing the analytical gradient outside the manifold. I overlooked it.

Because, I am modifying my point in the loop, but computing the analytical gradient before the assign.

[Nimrais]

Thanks for the fast response!

[kellertuer]

Ah indeed, you correct the point in the loop – then just move that fix (to bring p onto the manifold) before evaluating the analytical gradient and everything should look like in my post.

Great that we could fix it so fast – and thanks for using Manopt.jl