# Get started: Optimize!

This example illustrates how to set up and solve optimization problems and how to further get data from the algorithm using DebugOptions and RecordOptions. We will use the Riemannian mean and median as simple examples.

To start from the quite general case: A Solver is an algorithm that aims to solve

$\operatorname*{argmin}_{x∈\mathcal M} f(x)$

where $\mathcal M$ is a Manifold and $f:\mathcal M → ℝ$ is the cost function.

In Manopt.jl a Solver is an algorithm that requires a Problem p and Options o. While former contains static data, most prominently the manifold $\mathcal M$ (usually as p.M) and the cost function $f$ (usually as x->get_cost(p, x)), the latter contains dynamic data, i.e. things that usually change during the algorithm, are allowed to change, or specify the details of the algorithm to use. Together they form a plan. A plan uniquely determines the algorithm to use and provide all necessary information to run the algorithm.

## Example

A gradient plan consists of a GradientProblem with the fields M, cost function $f$ as well as gradient storing the gradient function corresponding to $f$. Accessing both functions can be done directly but should be encapsulated using get_cost(p,x) and get_gradient(p,x), where in both cases x is a point on the Manifold M. Second, the GradientDescentOptions specify that the algorithm to solve the GradientProblem will be the gradient descent algorithm. It requires an initial value o.x0, a StoppingCriterion o.stop, a Stepsize o.stepsize and a retraction o.retraction and it internally stores the last evaluation of the gradient at o.gradient for convenience. The only mandatory parameter is the initial value x0, though the defaults for both the stopping criterion (StopAfterIteration(100)) as well as the stepsize (ConstantStepsize(1.) are quite conservative, but are chosen to be as simple as possible.

With these two at hand, running the algorithm just requires to call x_opt = solve(p,o).

In the following two examples we will see, how to use a higher level interface that allows to more easily activate for example a debug output or record values during the iterations

## The given Dataset

using Manopt, Manifolds
using Random, Colors

For a persistent random set we use

n = 100
σ = π / 8
M = Sphere(2)
x = 1 / sqrt(2) * [1.0, 0.0, 1.0]
Random.seed!(42)
data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n]

and we define some colors from Paul Tol

black = RGBA{Float64}(colorant"#000000")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")
TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")

Then our data rendered using asymptote_export_S2_signals looks like

asymptote_export_S2_signals("startDataAndCenter.asy";
points = [ [x], data],
colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal]),
dot_size = 3.5, camera_position = (1.,.5,.5)
)
render_asymptote("startDataAndCenter.asy"; render = 2)

## Computing the Mean

To compute the mean on the manifold we use the characterization, that the Euclidean mean minimizes the sum of squared distances, and end up with the following cost function. Its minimizer is called Riemannian Center of Mass.

Note

There are more sophisticated methods tailored for the specific manifolds available in Manifolds.jl see mean.

F(M, y) = sum(1 / (2 * n) * distance.(Ref(M), Ref(y), data) .^ 2)
gradF(M, y) = sum(1 / n * grad_distance.(Ref(M), data, Ref(y)))

note that the grad_distance defaults to the case p=2, i.e. the gradient of the squared distance. For details on convergence of the gradient descent for this problem, see [Afsari, Tron, Vidal, 2013]

The easiest way to call the gradient descent is now to call gradient_descent

xMean = gradient_descent(M, F, gradF, data[1])

but in order to get more details, we further add the debug= options, which act as a decorator pattern using the DebugOptions and DebugActions. The latter store values if that's necessary, for example for the DebugChange that prints the change during the last iteration. The following debug prints

# i | x: | Last Change: | F(x):

as well as the reason why the algorithm stopped at the end. Here, the format shorthand and the [DebugFactory] are used, which returns a DebugGroup of DebugAction performed each iteration and the stop, respectively.

xMean = gradient_descent(
M,
F,
data[1];
debug=[:Iteration, " | ", :x, " | ", :Change, " | ", :Cost, "\n", :Stop],
)
Initial | x: [0.5737338264338113, -0.1728651513118652, 0.8005917410687816] |  | F(x): 0.22606088442202987
# 1 | x: [0.7823885991620455, 0.08556904440746377, 0.6168841208366815] | Last Change: 0.38188652881305873 | F(x): 0.14924728281088784
# 2 | x: [0.7917377109136624, 0.09720922045898238, 0.603076914311061] | Last Change: 0.020335997468986768 | F(x): 0.1490294899571422
# 3 | x: [0.792283271084327, 0.0977159913287459, 0.6022780117176446] | Last Change: 0.001092107177300155 | F(x): 0.14902886086292993
# 4 | x: [0.7923155575566623, 0.0977379000751398, 0.6022319819991105] | Last Change: 6.034190038620069e-5 | F(x): 0.14902885894006118
# 5 | x: [0.7923174650259398, 0.09773883592992506, 0.6022293205797368] | Last Change: 3.4056244498865836e-6 | F(x): 0.14902885893393145
# 6 | x: [0.7923175775065069, 0.09773887521445883, 0.6022291662199968] | Last Change: 1.9655981597711058e-7 | F(x): 0.14902885893391132
# 7 | x: [0.7923175841297168, 0.09773887682107818, 0.6022291572454814] | Last Change: 2.9802322387695312e-8 | F(x): 0.1490288589339113
# 8 | x: [0.792317584519285, 0.09773887688413527, 0.6022291567227154] | Last Change: 3.332000937312528e-8 | F(x): 0.14902885893391135
The algorithm reached approximately critical point after 8 iterations; the gradient norm (6.549992260414836e-10) is less than 1.0e-8.

A way to get better performance and for convex and coercive costs a guaranteed convergence is to switch the default ConstantStepsize(1.0) with a step size that performs better, for example the ArmijoLinesearch(). We can tweak the default values for the contractionFactor and the sufficientDecrease beyond constant step size which is already quite fast. This gives

xMean2 = gradient_descent(
M,
F,
data[1];
stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.5),
debug=[:Iteration, " | ", :x, " | ", :Change, " | ", :Cost, "\n", :Stop],
)
3-element Array{Float64,1}:
0.7923175844625908
0.09773887689188214
0.6022291567960479

which finishes in 5 steaps, just slightly better than the previous computation.

F(M, xMean) - F(M, xMean2)
3.3306690738754696e-16

Note that other optimization tasks may have other speedup opportunities.

For even more precision, we can further require a smaller gradient norm. This is done by changing the StoppingCriterion used, where several criteria can be combined using & and/or |. If we want to decrease the final gradient (from less that 1e-8) norm but keep the maximal number of iterations to be 200, we can run

xMean3 = gradient_descent(
M,
F,
data[1];
stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.5),
debug=[:Iteration, " | ", :x, " | ", :Change, " | ", :Cost, "\n", :Stop],
)
3-element Array{Float64,1}:
0.7923175845436088
0.09773887688651485
0.6022291566903287

which takes 10 iterations but gets a very small gradient, and not much is gained in the cost itself

F(M, xMean2) - F(M, xMean3)
1.3877787807814457e-16
asymptote_export_S2_signals("startDataCenterMean.asy";
points = [ [x], data, [xMean] ],
colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange]),
dot_size = 3.5, camera_position = (1.,.5,.5)
)
render_asymptote("startDataCenterMean.asy"; render = 2)

## Computing the Median

Note

There are more sophisticated methods tailored for the specific manifolds available in Manifolds.jl see median.

Similar to the mean you can also define the median as the minimizer of the distances, see for example [Bačák, 2014], but since this problem is not differentiable, we employ the Cyclic Proximal Point (CPP) algorithm, described in the same reference. We define

F2(M, y) = sum(1 / (2 * n) * distance.(Ref(M), Ref(y), data))
proxes = Function[(M, λ, y) -> prox_distance(M, λ / n, di, y, 1) for di in data]

where the Function is a helper for global scope to infer the correct type.

We then call the cyclic_proximal_point as

o = cyclic_proximal_point(
M,
F2,
proxes,
data[1];
debug=[:Iteration, " | ", :x, " | ", :Change, " | ", :Cost, "\n", 50, :Stop],
record=[:Iteration, :Change, :Cost],
return_options=true,
)
xMedian = get_solver_result(o)
values = get_record(o)
Initial | x: [0.5737338264338113, -0.1728651513118652, 0.8005917410687816] |  | F(x): 0.29868007939347
# 50 | x: [0.7897340334179458, 0.07829927185150307, 0.608431902918448] | Last Change: 0.08025247706840558 | F(x): 0.2454542345696772
# 100 | x: [0.7899680461384705, 0.0786029490954357, 0.6080888606722519] | Last Change: 0.0004925110933059937 | F(x): 0.24545405097297043
# 150 | x: [0.7900364131632669, 0.07868191382352603, 0.6079898208960564] | Last Change: 0.0001394128648418859 | F(x): 0.24545402791447266
# 200 | x: [0.7900676764403387, 0.07871541785855693, 0.6079448573970903] | Last Change: 6.2423013171665e-5 | F(x): 0.2454540216149762
# 250 | x: [0.7900851818616215, 0.07873311977053304, 0.607919814822658] | Last Change: 3.4399100660733445e-5 | F(x): 0.24545401919738047
The algorithm performed a step with a change (0.0) less than 1.0e-12.

where the differences to gradient_descent are as follows

• the third parameter is now an Array of proximal maps
• debug is reduces to only every 50th iteration
• we further activated a RecordAction using the record= optional parameter. These work very similar to those in debug, but they collect their data in an array. The high level interface then returns two variables; the values do contain an array of recorded datum per iteration. Here a Tuple containing the iteration, last change and cost respectively; see RecordGroup, RecordIteration, RecordChange, RecordCost as well as the RecordFactory for details. The values contains hence a tuple per iteration, that itself consists of (by order of specification) the iteration number, the last change and the cost function value.

values
254-element Array{Tuple{Int64,Float64,Float64},1}:
(1, 0.0, 0.24765709821769374)
(2, 0.03915137245266526, 0.24599998171726126)
(3, 0.014699094618037886, 0.24568197834841984)
(4, 0.007513039640263751, 0.24557305929190304)
(5, 0.004489348044805541, 0.24552491610823293)
(6, 0.002953601317333033, 0.24550007064482557)
(7, 0.0020748587400968685, 0.24548583414095682)
(8, 0.0015286431206780972, 0.24547703568610252)
(9, 0.0011677432888731137, 0.24547127727824178)
(10, 0.0009178169878496853, 0.2454673351043927)
⋮
(246, 3.18550442553044e-7, 0.2454540193270485)
(247, 2.9876735293512787e-7, 0.24545401929389457)
(248, 2.8585383434967613e-7, 0.24545401926124058)
(249, 2.682209014892586e-7, 0.2454540192290724)
(250, 2.2891598188485577e-7, 0.24545401919738047)
(251, 1.6858739404357632e-7, 0.24545401916616028)
(252, 1.1151007970493862e-7, 0.24545401913540696)
(253, 7.884953353001449e-8, 0.24545401910510958)
(254, 0.0, 0.24545401907525052)

The resulting median and mean for the data hence are

asymptote_export_S2_signals("startDataCenterMean.asy";
points = [ [x], data, [xMean], [xMedian] ],
colors=Dict(:points => [TolVibrantBlue, TolVibrantTeal, TolVibrantOrange, TolVibrantMagenta]),
dot_size = 3.5, camera_position = (1.,.5,.5)
)
render_asymptote("startDataCenterMedianAndMean.asy"; render = 2)`

## Literature

• [Bačák, 2014] Bačák, M: Computing Medians and Means in Hadamard Spaces., SIAM Journal on Optimization, Volume 24, Number 3, pp. 1542–1566, doi: 10.1137/140953393, arxiv: 1210.2145.
• [Afsari, Tron, Vidal, 2013] Afsari, B; Tron, R.; Vidal, R.: On the Convergence of Gradient Descent for Finding the Riemannian Center of Mass, SIAM Journal on Control and Optimization, Volume 51, Issue 3, pp. 2230–2260. doi: 10.1137/12086282X, arxiv: 1201.0925