Getting Started with STOP

Introduction

The STOP toolbox is designed for optimization problems on the Stiefel manifold, which could be expressed as

where $I_p$$I_p$ refers to the $p$$p$-th order identity matrix, $X$$X$ is a matrix with $n$$n$ rows and $p$$p$ columns. The feasible set of this optimization problem

can be regarded as a Riemannian manifold in ${\mathbb{R}}^{n\times p}$$\mathbb{R}^{n\times p}$, and we also call it Stiefel manifold.

This toolbox is an integration of several state-of-the-art solvers for optimization problems on the Stiefel manifold. STOP is a prototyping toolbox, designed by capturing the special structures of problems on the Stiefel manifold. For basic use, one only needs to describe the cost function and its derivatives and call the build-in solvers. Basically, all the options in STOP toolbox have default values, and one can execute a solver without any specified option. This toolbox is under maintenance and working in progress. Any feedbacks and contributions are welcome!

Install STOP

Install

The STOP toolbox can be downloaded from Gitee. The current version of STOP is STOPv0.1 and was packaged on 2021/10/12. The STOP toolbox supports two different ways of installation. One could directly copy the solvers to the working folder for temporary execution. Besides, One can also install STOP to Matlab by the following steps.

1. Extract the STOP.zip file to the installation folder, say /my_install_folder/.
2. Go to/my_install_folder/ folder and execute /my_install_folder/STOP/install_stop.m in Matlab. If you wish to save the path for future usage, you can use the savepath function to save the path to STOP solvers.
3. Enjoy!

Check

You can execute the script /my_install_folder/STOP/check_install/checkinstall.m to check whether STOP is properly installed. This script provides a overall check for all the solvers in STOP toolbox. If there are no errors, you are done! Otherwise, please feel free to contact us.

Example

Problem formulation

In this section, we consider the following nonlinear eigenvalue problem

where $\rho =\mathrm{Diag}(X{X}^{\mathrm{\top }})$$\rho = \mathrm{Diag}(XX^\top)$, and $L^{†}$$L^{\dagger}$ denotes the pseudo-inverse of the positive definite matrix $L$$L$, i.e. $L^{†}L{L}^{†}=L^{†}$$L^{\dagger}LL^{\dagger} = L^{\dagger}$, $L{L}^{†}L=L$$LL^{\dagger}L = L$. Here we uses $\mathrm{Diag}(M)$$\mathrm{Diag}(M)$ to denote the vector that is composed of diagonal entries of the square matrix $M$$M$, while $\mathrm{diag}(v)$$\mathrm{diag}(v)$ refers to a diagonal matrix with $v$$v$ being its diagonal entries. Then the cost function and its Euclidean gradient can be expressed as

In this example, we choose $L$$L$ as a tri-diagonal matrix generated by L = gallery('tridiag',n,-1,2,-1). Noting that $L$$L$ is full-rank, then we can conclude that $L^{†}=L^{-1}$$L^{\dagger} = L^{-1}$ in this case. We solve this simple optimization problem using solvers in STOP to illustrate the most basic usage of the STOP toolbox.

An Example Code

Solving this optimization problem using solvers in STOP requires little Matlab codes:

32
1
% Problem setting
2
n = 1000;  p = 10; 
3
L = gallery('tridiag',n,-1,2,-1); 
4
alpha = 1; 
5
[Ll,Lu] = lu(L);
6

7

8
% Set options
9
opts = [];  % no predefined parameter
10
opts.X = [];
11
opts.dim = [n,p]; %specify dimension
12

13
% Solve
14
t1 = tic;
15
[Out] = SLPG(@myfun,opts,L, Ll,Lu,alpha);
16
t=toc(t1)
17

18
% Plot
19
pic = semilogy(1:Out.iter, Out.kkts,'.-', 1:Out.iter, Out.feas,'*-');
20
xlabel('Number of iterations')
21
legend('Substationarity','Feasibility')
22

23

24
% Define objective function
25
function [funX, F] = myfun(X, L, Ll,Lu,alpha)
26
    LX = L*X;
27
    rhoX = sum(X.^2, 2); 
28
    tempa = Lu\(Ll\rhoX); 
29
    tempa = alpha*tempa;
30
    F = LX + bsxfun(@times,tempa,X);
31
    funX = 0.5*sum(sum(X.*(LX))) + 1/4*(rhoX'*tempa); 
32
end

Let us review the codes step by step. First, we generate the data (matrix $L$$L$) for the optimization problem by the following codes.

4
1
n = 1000;  p = 10; 
2
L = gallery('tridiag',n,-1,2,-1); 
3
alpha = 1; 
4
[Ll,Lu] = lu(L);

Here instead of directly computing $L^{-1}\rho$$L^{-1}\rho$ in each step, we choose to factorize $L=Ll\ast Lu$$L = Ll*Lu$ by LU factorization at the beginning stage of the codes. Therefore, in each step we could compute Lu\(Ll\rho) to compute $L^{-1}\rho$$L^{-1}\rho$, which could be more efficient.

Then we specify the cost function and its gradient in the following function.

9
1
% Define objective function
2
function [funX, F] = myfun(X, L, Ll,Lu,alpha)
3
    LX = L*X;
4
    rhoX = sum(X.^2, 2); 
5
    tempa = Lu\(Ll\rhoX); 
6
    tempa = alpha*tempa;
7
    F = LX + bsxfun(@times,tempa,X);
8
    funX = 0.5*sum(sum(X.*(LX))) + 1/4*(rhoX'*tempa); 
9
end

Currently, in STOP toolbox, we require the function return the function value and its gradient simultaneously. Usually, computing the function value and gradient simultaneously is much faster than compute them separately, even when cache techniques are involved. To achieve a better performance, we strongly suggest to compute the function value and gradient in a single function.

Now, with cost function and the data, we can set options for the solver. The initial point or dimensions should be specified in opts,

3
1
opts = [];  % no predefined parameter
2
opts.X = [];
3
opts.dim = [n,p]; %specify dimension

Then we call a solver, for instance, SLPG to solve the nonlinear eigenvalue problem.

3
1
t1 = tic;
2
[Out] = SLPG(@myfun,opts,L,Ll,Lu,alpha);
3
t = toc(t1)

Here the cost function and its derivatives are called by function handles.

Finally, we access the contents of the output structure Out and display the convergence plot of our solver.

3
1
pic = semilogy(1:Out.iter, Out.kkts,'.-', 1:Out.iter, Out.feas,'*-');
2
xlabel('Number of iterations')
3
legend('Substationarity','Feasibility')

Solvers

Solvers, or optimization algorithms, are functions in STOP.

In principle, all solvers admit the basic call format as follows.

1
1
[out] = solver(fun,options);

And we could only specify options.dim in options.

The output structure out contains the final solution to the optimization problem that can be accessed by out.X. Besides, the output structure out may contain the information to describe the optimality of the final solution (including the F-norm of Riemannian gradient, function value, and feasibility at the last iterate).

As Stiefel manifold is a nonconvex set in ${\mathbb{R}}^{n\times p}$$\mathbb{R}^{n\times p}$, the optimization problems on Stiefel manifold is usually nonconvex and the final solution out.X is not guaranteed to be a global minimizer.

Furthermore, all solvers also admit a more complete call format.

1
1
[out] = solver(fun,options,varargin);

Here the output structure out may contain information collected in the execution of the solver, including the time, feasibility, substationarity.

The input options options enables a customized parameters for the provided solver, including the parameters that affects the performance of the solver, the stopping criteria, and the information to be displayed in the execution of the solver.

Back to Homepage