Quasi-Newton and Conjugate Gradient on Stiefel/Unitary/Orthogonal Manifolds

[BenjaminWoeckinger]

Hello!
My aim is to perform quasi-newton and conjugate gradient algorithms(preferably Polak Ribiere) on the Unitary and Orthogonal manifolds. Sadly, for Orthogonal and Unitary Matrices, parallel transport is apparently not implemented, although it is set as default vector transport method (in theory, there would be a closed form for the unitary/orthogonal manifold though).

I tested some configurations on the Brockett optimization problem (a diagonalization problem which is also featured in your github tests (test_quasi_Newton.jl)) and ran into some problems.
I didn't get my code running for Unitary Matrices, because it doesn't get along with complex numbers (errors are: "no method matching isless(::float64,::complexF64)" or "log!(::AbstractManifold, ::Matrix{ComplexF64}, ::Matrix{ComplexF64}, ::Matrix{ComplexF64}) is ambiguous"). I'm also not really familiar with julia, so I'm having troubles finding the corresponding code. This problem also occurs with the complex Stiefel manifold. I'm not sure, if that's a fixable problem or just not supported at the moment.

In general, I don't really get the CG algorithm going. Polak-Ribiere doesn't work at all, Fletcher Reeves works sometimes, depending on the initial conditions and size of the problem.

At times, I also have problems with the log function (leading to "matrix contains Infs/NaNs" errors), but this is not always reproducible. For the real Stiefel manifold, some log/exp tests fail when using test_manifold for larger matrices! (10x10 for example) Maybe the error tolerance has to be set larger here, I did not go too much into detail.

Most importantly, which retractions, vector transports and stepsize algorithms would you recommend for these kind of problems? Ideally, I would like to work on the Unitary Matrices manifold with the exponential map, parallel transport and a polynomial linesearch, but I would need to include that to the module I guess. Do you think it would bring a performance boost, or are other methods similar in accuracy and performance?

This is the code I am performing my tests with at the moment (it is less complicated than it looks at first glance):

using Test
using LinearAlgebra
using Manifolds
using ManifoldsBase
using ManifoldDiff
using Random
using Manopt

################################ set parameters
        algo = "BFGS"                   #set algorithm (BFGS or CG)
        BFGS_memory  = 10               #set memory size of BFGS algorithm
        manifold = "Orthogonal"         #choose between Stiefel, Unitary and Orthogonal
        complex_calc = false            #set if manifold shall be complex (automatically set for Unitary and orthogonal)
        dim = 5 
        random_start = true             #if random_start=false, dim must be 4 and complex_calc must be false here

############################### define manifold
        if manifold == "Stiefel"
                if complex_calc
                        M = Stiefel(dim,dim,ℂ)
                else
                        M = Stiefel(dim,dim)
                end
        elseif manifold == "Unitary"
                M = UnitaryMatrices(dim)
                complex_calc = true
        elseif manifold == "Orthogonal"
                M = OrthogonalMatrices(dim)
                complex_calc = false
        end

############################## set Matrix for Brocket optimization problem
        N = Diagonal{Float64}(1:1:dim)

############################### define starting conditions (hermitian and unitary)
        if random_start
                if complex_calc
                        H = Hermitian(rand(dim,dim).-0.5+(rand(dim,dim).-0.5)*im) #creates a random complex Hermitian
(+-0.5)
                        U = rand(Stiefel(dim,dim,ℂ))                            #creates a random unitary matrix
                else
                        H = Hermitian(rand(dim,dim).-0.5)                       #creates a ramdom Hermitian (+-0.5)
                        U = rand(Stiefel(dim,dim))                              #creates a random orthogonal matrix
                end
        else
                H = Hermitian([2.0 1.0 0.0 3.0; 1.0 3.0 4.0 5.0; 0.0 4.0 3.0 2.0; 3.0 5.0 2.0 6.0]) #dim has to be 4!
                U = Matrix{Float64}(I, dim, dim)                                #only real for now!
        end
        
################################ define costfunction
        f(M,x) = real(tr(x' * H * x * N))

################################# define gradient (left translation of gradient for Unitary/orthogonal)
        function grad_f(M,z)
                Γ = H * z * N 
                if ((manifold == "Unitary") || (manifold == "Orthogonal"))
                        return (z' * Γ -  Γ' * z)
                elseif manifold == "Stiefel" 
                        return (Γ - z * Γ' * z)
                end
        end
################################# perform optimization
        if algo == "CG"
                U = conjugate_gradient_descent(M, f, grad_f, U; 
#                       stepsize = ArmijoLinesearch(M; retraction_method = QRRetraction()), 
                        stepsize = ArmijoLinesearch(M),
                        retraction_method = QRRetraction(),
                        vector_transport_method = ProjectionTransport(),
#                       coefficient = PolakRibiereCoefficient(ProjectionTransport()), #doesn't work properly
#                       coefficient = SteepestDirectionUpdateRule(),   #works kind of (bad convergence obv)
                        coefficient = FletcherReevesCoefficient(),   #works sometimes (better for small dim)       
                        stopping_criterion = StopAfterIteration(5000) |StopWhenGradientNormLess(1e-10),
            
                        debug=[:Iteration,(:Change, "|ΔU|: %1.9e |"),(:GradientNorm, " |G|: %1.9e |"),
                                (:Cost, " F(x): %1.11e | "), "\n", :Stop]
                )

        elseif algo == "BFGS"
                U = quasi_Newton(M, f, grad_f, U;
                        memory_size = BFGS_memory,
                        stepsize = WolfePowellBinaryLinesearch(
                                                             M;
#                                                            retraction_method = QRRetraction(),
                                                             vector_transport_method = ProjectionTransport(),
                                                             ),
                        vector_transport_method = ProjectionTransport(),
#                       retraction_method = QRRetraction(),
                        cautious_update = true,
                        stopping_criterion = StopWhenGradientNormLess(1e-10),
                        debug = [:Iteration,(:Change, "|ΔU|: %1.9e |"),(:GradientNorm, " |G|: %1.9e |"),
                                (:Cost, " F(x): %1.11e | "), "\n", :Stop]
                )
        end

Maybe you can help me with this, thank you in advance,
Benjamin

[kellertuer]

Hi,
thanks for using Manopt.jl.

I can check your method closer the next few days, but with so many lines of code, that really will take a bit of time, reading and following up on so much code – even if you say it is simpler than it looks – at first glance it is more like a wall-of-code.
A Minimal (non-=Working example would help more – and a precise error message (including trace.

A first few comments

Sadly, for Orthogonal and Unitary Matrices, parallel transport is apparently not implemented, although it is set as default vector transport method (in theory, there would be a closed form for the unitary/orthogonal manifold though).

This is more a matter of Manifolds.jl.
We do our best (11k lines of code, even more in docs) to implement as much as we are aware of on manifolds.
If no retraction is implemented and no default is set, the global default is ParallelTransport() (in the idea of: by default we hope for the best).
The reason that they are (maybe) known but not implemented might simply have the reason, that no one found the time (or reference) yet to implement (document and test) these. Please open an issue than, for best of cases with the closed form solution linked, then we can see what we can do.
On the other hand, this is currently done by two people in their free time, so also open issues might take a bit every now and then. But we are always interested in more ideas.

Compared to that, the work on Manopt.jl is even only one person mainly developing that – me.

In general, I don't really get the CG algorithm going. Polak-Ribiere doesn't work at all, Fletcher Reeves works sometimes, depending on the initial conditions and size of the problem.

This is a bit hard to help, if the problem is a bit vague in “could not get it to work”

For the real Stiefel manifold, some log/exp tests fail when using test_manifold for larger matrices! (10x10 for example) Maybe the error tolerance has to be set larger here, I did not go too much into detail.

Yes that is probably a tolerance issue

Most importantly, which retractions, vector transports and stepsize algorithms would you recommend for these kind of problems? Ideally, I would like to work on the Unitary Matrices manifold with the exponential map, parallel transport and a polynomial linesearch, but I would need to include that to the module I guess. Do you think it would bring a performance boost, or are other methods similar in accuracy and performance?

The retraction / vector transport depends a lot on the problem you have, the manifold size you habe and such. But for unitary both exp ( https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/generalunitary.html#Base.exp-Tuple{Manifolds.GeneralUnitaryMatrices,%20Any,%20Any} ) and log (https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/generalunitary.html#Base.log-Tuple{Manifolds.GeneralUnitaryMatrices,%20Any,%20Any}, maybe with a typo in the documentation it seems) are available, if they are, maybe try them first over retractions.
Parallel transport might be a nice addition to Manifolds.jl as far as I can see.
If either of these struggles with accuracy on large scale cases – it might be worth checking the retractions and vector transports as well.

The step size depends a lot on the problem at hand and the algorithm you want to use. For example „a polynomial linesearch“ is currently probably not available but feel free to contribute that. Since I have note yet worked with that I can not say anything about possible performance boost (of things that I did not yet use).

[BenjaminWoeckinger]

Hi, thanks for your quick response and your effort in general!

I can check your method closer the next few days, but with so many lines of code, that really will take a bit of time, reading and following up on so much code – even if you say it is simpler than it looks – at first glance it is more like a wall-of-code.
A Minimal (non-=Working example would help more – and a precise error message (including trace.

That would be amazing, thank you! Sorry for the mess, it's the options for all three manifolds (and random initialization) that blow up the code a bit. Here is a shorter version thats more understandable, I hope. For complete error messages, I'll come back to you, I will be on holidays next week and need to trace down the issues more systematically!

The problem is basically the same as "Brocket" in "test/solvers/test_quasi_Newton.jl", but simpler because I have only square matrices. (I tried three manifolds, unitary, orthogonal and Stiefel as their generalization though). A Hermitian matrix H shall be diagonalized by optimizing the Brockett functional, using a diagonal matrix N with distinct entries.

In the end, the important things are the two blocks U1 and U2, where the optimization takes place.
For the quasi-Newton method, I used mostly the same keywords, so I thought maybe similar problems emerged and these parameters (and not the default values) were chosen for specific reasons that I'm not aware of.

dim=10

M = Stiefel(dim,dim)

#M = Stiefel(dim,dim,ℂ)           
#M = UnitaryMatrices(dim)
#M = OrthogonalMatrices(dim)

############################## set Matrix for Brocket optimization problem
N = Diagonal{Float64}(1:1:dim)

############################## set starting point U and Hermitian H
U = rand(M)                         
H = Hermitian(rand(dim,dim).-0.5)

################################ define costfunction
f(M,x) = real(tr(x' * H * x * N))

################################# define gradient
function grad_f(M,x)
   Γ = H * x * N 
   return (Γ - x * Γ' * x)
end

################################# perform optimization
U1 = conjugate_gradient_descent(M, f, grad_f, U; 
                        stepsize = ArmijoLinesearch(M),
                        retraction_method = QRRetraction(),
                        vector_transport_method = ProjectionTransport(),
#                       coefficient = PolakRibiereCoefficient(ProjectionTransport()), #doesn't work properly
#                       coefficient = SteepestDirectionUpdateRule(),   #works kind of (bad convergence obv)
                        coefficient = FletcherReevesCoefficient(),   #works sometimes (better for small dim)       
                        stopping_criterion = StopAfterIteration(5000) |StopWhenGradientNormLess(1e-10),
            
                        debug=[:Iteration,(:Change, "|ΔU|: %1.9e |"),(:GradientNorm, " |G|: %1.9e |"),
                                (:Cost, " F(x): %1.11e | "), "\n", :Stop]
                )

U2 = quasi_Newton(M, f, grad_f, U;
                        memory_size = BFGS_memory,
                        stepsize = WolfePowellBinaryLinesearch(
                                                             M;
#                                                           retraction_method = QRRetraction(),
		                                             vector_transport_method = ProjectionTransport(),
                                                             ),
                        vector_transport_method = ProjectionTransport(),
#                       retraction_method = QRRetraction(),
                        cautious_update = true,
                        stopping_criterion = StopWhenGradientNormLess(1e-10),
                        debug = [:Iteration,(:Change, "|ΔU|: %1.9e |"),(:GradientNorm, " |G|: %1.9e |"),
                                (:Cost, " F(x): %1.11e | "), "\n", :Stop]

This is more a matter of Manifolds.jl.
We do our best (11k lines of code, even more in docs) to implement as much as we are aware of on manifolds.
If no retraction is implemented and no default is set, the global default is ParallelTransport() (in the idea of: by default we hope for the best).
The reason that they are (maybe) known but not implemented might simply have the reason, that no one found the time (or reference) yet to implement (document and test) these. Please open an issue than, for best of cases with the closed form solution linked, then we can see what we can do.
On the other hand, this is currently done by two people in their free time, so also open issues might take a bit every now and then. But we are always interested in more ideas.
Compared to that, the work on Manopt.jl is even only one person mainly developing that – me.

Wow thats really cool! I understand, I just wanted to point it out, I thought maybe it's an error that nobody is aware of.

This is a bit hard to help, if the problem is a bit vague in “could not get it to work”

That's true, I was a bit lost here and just tried some coefficients with several retractions and vector transport methods, but most did not converge or only very slowly.
The cost function we want to optimize eventually, is a polynomial in the elements u_ij of a unitary matrix U.
Do you have an educated guess what might work best for this manifold/problem that can be done in manopt? We already have a CG implementation with the Polak Ribiere Coefficient that works quite well for comparison (following this paper, just if you are interested: https://www.cs.ox.ac.uk/files/6671/AbrEriKoi09SP__CG_for_optimization_under_unitary_matrix_constraints_appeared.pdf) and I was curious if I could reproduce/compare it with your module.

In the end, my goal is, to get a (limited memory) BFGS algorithm that can outperform our existing CG algorithm. I did program a L-BFGS already (following https://www.math.fsu.edu/~aluffi/archive/paper488.pdf), but it didn't perform as expected. In particular, more memory lead to worse performance, which is why I wanted to check for consistency and maybe also get a comparison to the "full" BFGS.

Thanks a lot for your help!

[kellertuer]

Or to make the distinction more precise – Manifolds.jl does have an AbstractLieGroup concept, but that is relatively new and not yet in `ManifoldsBase.jl (the interface Manopt is based on) – so in Manopt it is not possible yet to write things purely for Lie groups.

BUt at least all this explains why for now the code did not work well – you took the manifold not the Lie group, and especially tangent vectors are tremendously different (especially in how they are represented/stored) for these two things.

[kellertuer]

I just found a nice fix. UnitaryMatrices and hence Unitary use the Lie algebra, so it was easy to adapt rand and project to this setting. On both you can soon (see JuliaManifolds/Manifolds.jl#706) run your experiment and the gradient you provide is correct. on Stiefel however, we have the representation I described here, so there you have to take your gradient times the base point you are at to make it correct.

[kellertuer]

This one is finished and merged, I hope that help :)

[BenjaminWoeckinger]

Also things “just” on Lie groups are currently a bit complicated to realise in Manopt I fear – and I would not have time to implement such an advanced algorithm before – uff – July? This current time is called “semester”.

Ok, I understand, then it's most likely not feasible for us to implement on our own. Would it help if we provided the code for this? We have it written in julia.

Thanks for the fixes! (I'll answer to the Lie stuff below)

[kellertuer]

For now with the new randoms it should hopefully work (since we removed the error message above).

[mateuszbaran]

Hi!

I didn't get my code running for Unitary Matrices, because it doesn't get along with complex numbers (errors are: "no method matching isless(::float64,::complexF64)" or "log!(::AbstractManifold, ::Matrix{ComplexF64}, ::Matrix{ComplexF64}, ::Matrix{ComplexF64}) is ambiguous").

@kellertuer I think one problem here is that this line:

Manopt.jl/src/solvers/quasi_Newton.jl

Line 658 in ebd4833

if sk_normsq != 0 && (inner(M, p, st.sk, st.yk) / sk_normsq) >= bound

should take only the real part of the inner product.

How did you get the log! issue? There might be an easy fix for that but it's hard to tell if I don't know how to reproduce it.

Most importantly, which retractions, vector transports and stepsize algorithms would you recommend for these kind of problems? Ideally, I would like to work on the Unitary Matrices manifold with the exponential map, parallel transport and a polynomial linesearch, but I would need to include that to the module I guess. Do you think it would bring a performance boost, or are other methods similar in accuracy and performance?

I would recommend trying the new Hager-Zhang linesearch: https://github.com/mateuszbaran/ImprovedHagerZhangLinesearch.jl . It performs much better than WolfePowell in my experiments. There are still a few things I need to tweak there and add an installation instruction but I will try to do it in the next couple of days.

Regarding retractions and VTs, it's very much problem-dependent. I'm currently working on an automatic tuner that chooses them: https://github.com/JuliaManifolds/Manopt.jl/pull/339/files#diff-ee1b3e2c8f8de306b636f704ead86f0db8063ab3af536fdf14239c7d85c9ca52 . If you are interested, I can help you try that on your problem. It uses a Bayesian hyperparameter optimization with early stopping, using Python's Optuna library (Julia doesn't have anything close in functionality).

In the end, my goal is, to get a (limited memory) BFGS algorithm that can outperform our existing CG algorithm. I did program a L-BFGS already (following https://www.math.fsu.edu/~aluffi/archive/paper488.pdf), but it didn't perform as expected. In particular, more memory lead to worse performance, which is why I wanted to check for consistency and maybe also get a comparison to the "full" BFGS.

L-BFGS is quite hard to get right but I would expect it to work better than CG. It's definitely a good idea to investigate.

[mateuszbaran]

I'm not sure about gradient stuff but Hager-Zhang is already implemented here: https://github.com/mateuszbaran/ImprovedHagerZhangLinesearch.jl/ . You can just

(@v1.10) pkg> dev https://github.com/mateuszbaran/ImprovedHagerZhangLinesearch.jl/

In your Julia environment and then pass HagerZhangLinesearch(M) as the line search (or stepsize) for optimization procedure.

[kellertuer]

That is already great news!

For the check gradient stuff, it doesn't pass the test consistently, even at quite high tolerance values (0.1 for example)
Still, oftentimes it passes at very small tolerances (1e-12) which speaks against a simple numerical accuracy issue.

Teh gradient check is most probably not the cause here. That is adapted code from the much more mature Matlab package. If you load Plots.jl you can even activate a plot on which scales the gradient is wrong, but with so large errors, I would assume, that your gradient is only sometimes correct.

What I also noticed, BFGS sometimes just freezes

This is most probably also the reason for BFGS (great that it works now!) freezes.
my interpretation that it happens more often in larger systems is that, there the probability that you “hit” a point where your gradient is wrong is maybe also higher.
With a correct gradient BFGS should not do that.

Stiefel seems to be a bit more stable than Orthogonal.

I have not yet had the chance to compare the two, so I can not comment much on that.

[kellertuer]

I took a bit of time and used your code above but for the example above, all works fine.

I did generate a Pluto notebook, but one can only attach txt files, so one has to remove the .txt from this and use that as a start
check_BW.jl.txt

Could we maybe try to generate a case that fails? This one for example does the following gradient check:

which indicates that for very small step sizes, it might get a bit imprecise, for best of cases the slope should be 2, here it is not completely, but the main reason is the very small steps on the left.

The code runs but takes a while since the precision 10e-11 is a bit high for Stiefel to be fast (I think).

If you could provide a similar example that fails (either with gradient or with conversion) I can surely take a look :)

[BenjaminWoeckinger]

In your Julia environment and then pass HagerZhangLinesearch(M) as the line search (or stepsize) for optimization procedure.

Thanks, this apparently solved the issue with the algorithm freezing! (Before I was using WolfePowellBinaryLinesearch)

[BenjaminWoeckinger]

How does this Pluto notebook thing work exactly? Can we work on the same file simultaneously? I only made a local copy for now:
check_BW_2.jl.txt
This fails most of the time for example (apart from not being able to install the Hager Zhang package, I don't know if it's possible there)

[kellertuer]

Nice! If you want to be on the safe side, check the gradient (with the plot for example), but Yes WolfePowellBinary might get stuck, I experienced that only once until now, but that is also hard to track down, the algorithm behind is not so easy to debug.

Great that it seems to work now!

If you end up with a nice numerical experiment maybe consider adding it to the ManoptExamples.jl package :)

[BenjaminWoeckinger]

Thanks for looking into it!

[BenjaminWoeckinger]

What does the check_gradient plot show in particular? is it the error plotted over the gradient norm?

[kellertuer]

I just saw that in the beginning you wrote

Sadly, for Orthogonal and Unitary Matrices, parallel transport is apparently not implemented, although it is set as default vector transport method (in theory, there would be a closed form for the unitary/orthogonal manifold though).

Are these still missing? If so, could you maybe point to the closed form you spoke about so we can add them to Manifolds.jl :) ?

[BenjaminWoeckinger]

Yes it's still missing, but is it needed if I use the group and not the manifold? After all, the vectors aren't really transported there (I still have to choose a different transport method there though to not get an error)
Sure, a closed form is here: https://doi.org/10.48550/arXiv.physics/9806030 at page 311
or here https://diginole.lib.fsu.edu/islandora/object/fsu:185265/datastream/PDF/view in slightly different notation, but I think their source is the Edelman paper above

[kellertuer]

Ok, we could check at some point to add that – just that I will not have much time before end of May (end of corrections for the exams).

The reason it does return ParallelTransport, is that that is default in the sense that: This way all that do implement it, do not have to set the default_vector_transport_method to indicate that they did implement it

[BenjaminWoeckinger]

Projection transport is working, maybe it would be convenient for others to just put that as default for the moment?

[kellertuer]

Oh, sometimes I miss a bit of overview. Indeed, then we should put that as a default. Do you want to do a PR?

[mateuszbaran]

I can handle a PR. Currently we have issues with CI so merging it may take some time.