-
Notifications
You must be signed in to change notification settings - Fork 41
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
trust_regions loses normalization on Sphere #306
Comments
Hi,
julia> 1-1.0000000000000002
-2.220446049250313e-16 so it might even just be the dimension of the sphere causing an error in the computation of the norm, since this is of order of machine precision. BackgroundThere is actually 2 things interacting here
The fix could be to always project the sub solver result (back) one the tangent space. The disadvantage is, that that might be costly on some manifolds or sometimes even not be implemented or defined. But we have therfore the optional Updated ExampleHm I could not directly run your example since with your packages (at least on my Julia) using Random, Manifolds, Manopt, Zygote, SparseArrays
# we set up a simple function to keep the optimizer busy.
n = 100
Random.seed!(0)
A = sprandn(n, n, 0.2); A = A + A'
M = Manifolds.Sphere(n-1, ℂ)
x0 = project(Manifolds.Sphere(size(A,1) - 1, ℂ), randn(Complex{eltype(A)}, size(A, 1)))
# this has no minimum since A is not posdef, but we don't care.
f(M, x) = real(x'*A*x)
g(M, x) = project(M, x, gradient(x->f(M,x), x)[1]) # you could also check out riemannian_gradient
x_old = trust_regions(M, f, g, x0)
x = trust_regions(M, f, g, x0; (project!)=project!) then we have julia> norm(x)
1.0000000000000027 so that is in the area of and we even see
Not only is the new cost lower (ok the old point is also not in the domain of the “shortened version of Rayleigh” (no division by the norm); but also that the new one does have a small gradient. Slightly relatedI started reworking trust regions to finally allow other sub solvers – but that is still a bit of work (and needs the next breaking version of ManifoldsBase.jl, 0.15) see #294, since we discussed that earlier. So this might still take some time, but will be available some time. One more thingHere, with projection the algorithm is also much faster (ok because it stops with a small gradient not due to max iterations ;)) julia> @btime x = trust_regions($M, $f, $g, $x0; (project!)=$project!);
5.454 ms (12737 allocations: 41.69 MiB)
julia> @btime x_old = trust_regions($M, $f, $g, $x0);
102.932 ms (219532 allocations: 636.69 MiB) |
Thanks for the very detailed answer! I missed that option in the docs, sorry. I confirm also my 'true' function of interest looks a bit like a Rayleigh quotient. And I confirm I forgot to include |
Great to here that that solved it :) I agree that by default not having the project (and hence these cases) might be a point of discussion, but it is a choice between speed and availability on more manifolds (like fixed rank) versus stability/inaccuracy for a few others. If you feel we could improve the docs here, let me know. |
The
and this problem seems to be specific to the I could submit a PR myself if you wish, but of course I am far from being an expert since I did not even know about this issue until last week. :) Many thanks for your work on this package! |
I am rework ing trust regions here #294 anyways, But I will add such a comment for sure, no problem :) Thanks for the feedback! |
After running sufficiently many iterations of
trust_regions
on the complex sphere, the result no longer has norm 1. Here is a minimal example:Returns
On more complicated examples, the norm gets even farther away from 1.
The text was updated successfully, but these errors were encountered: