# PySTOP
## Introduction
The STOP toolbox is designed for **optimization problems on the Stiefel manifold**, which could be expressed as
$$
\begin{aligned}
\min_{X \in \mathbb{R}^{n\times p}} ~ &f(X)\\
\text{s. t.}~& X^\top X = I_p,
\end{aligned}
$$
where $I_p$ refers to the $p$-th order identity matrix, $X$ is a matrix with $n$ rows and $p$ columns. The feasible set of this optimization problem
$$
\mathcal{S}_{n,p} := \left\{X \in \mathbb{R}^{n\times p}: X^\top X = I_p \right\},
$$
can be regarded as a Riemannian manifold in $\mathbb{R}^{n\times p}$, and we also call it as **Stiefel manifold**.
This document describes the python version of the STOP package (PySTOP). Currently, PySTOP only involves the SLPG algorithm, which could handle smooth, $\ell_1$-norm regularized and $\ell_{2,1}$-nom regularized optimization problems on the Stiefel manifold. The package is
## Installation
The source code of PySTOP package can be found from [the website](https://stmopt.gitee.io/). Besides, it supports direct installation from `pip`:
```shell
pip install pystop
```
## Example
### Problem formulation
In this section, we consider the following nonlinear eigenvalue problem
$$
\min_{X \in \mathcal{S}_{n, p}} ~ \frac{1}{2}\mathrm{tr}(X^\top L X) + \frac{\alpha}{4} \rho^\top L^{\dagger} \rho,
$$
where $\rho = \mathrm{Diag}(XX^\top)$, and $L^{\dagger}$ denotes the pseudo-inverse of the positive definite matrix $L$, i.e. $L^{\dagger}LL^{\dagger} = L^{\dagger}$, $LL^{\dagger}L = L$. Here we uses $\mathrm{Diag}(M)$ to denote the vector that is composed of diagonal entries of the square matrix $M$, while $\mathrm{diag}(v)$ refers to a diagonal matrix with $v$ to be its diagonal entries. Then the cost function and its **Euclidean gradient** can be expressed as
$$
\begin{aligned}
f(X) ={}& \frac{1}{2}\mathrm{tr}(X^\top L X) + \frac{\alpha}{4} \rho^\top L^{\dagger} \rho,\\
\nabla f(X) ={}& LX + \alpha \mathrm{diag}(L^{\dagger}\rho)X.
\end{aligned}
$$
In this example, we choose $L$ as a tri-diagonal matrix generated by `L = gallery('tridiag',n,-1,2,-1)`. Noting that $L$ is full-rank, then we can conclude that $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. For additional theory, readers are recommended to refer the papers in the about page.
```python
# Import packages
import numpy as np
import scipy as sp
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
# Import manifolds and solvers
from pystop.manifold import Stiefel
from pystop.solver import SLPG_smooth
# Set parameters
n = 1000
p = 10
alpha = 1
M = Stiefel(n,p)
# Defining objective function
L = diags(np.array([-1, 2, -1]), np.array([1, 0, -1]), shape = (n,n)).tocsc()
def obj_fun(X):
LX = L@X
rho = np.sum(X * X, 1)
Lrho = spsolve(L, rho)
fval = 0.5*np.sum(X* LX) + (alpha /4) * np.sum(rho * Lrho)
grad = LX + alpha * Lrho[: ,np.newaxis] * X
return fval, grad
# Execute the solver
X, out_dict = SLPG_smooth(obj_fun, M)
```
Let us review the code step by step. First, we specify the dimension of the problem and specify the Stiefel manifold. In `pySTOP` package, we need to specify the dimension of the Stiefel manifold before executing the solver. The Stiefel manifold should be specified as the *STOP manifold class*, for example,
```python
# Set parameters
n = 1000
p = 10
alpha = 1
# Specify the Stiefel manifold
M = Stiefel(n,p)
```
Here `pystop.manifold.stiefel` is a build-in function to specify the Stiefel manifold and hence provides essential tools for the algorithm.
Then we generate the data (matrix $L$) for the optimization problem by the following code,
~~~python
L = diags(np.array([-1, 2, -1]), np.array([1, 0, -1]), shape = (n,n)).tocsc()
~~~
Here we utilize `SciPy.sparse` to create a sparse representation of $L$ . Therefore, in each step we could use the `scipy.sparse.linalg.spsolve` function to compute .
Then we specify the cost function and its gradient in the following function
~~~python
# Defin objective function
def obj_fun(X):
LX = L@X
rho = np.sum(X * X, 1)
Lrho = spsolve(L, rho)
fval = 0.5*np.sum(X* LX) + (alpha /4) * np.sum(rho * Lrho)
grad = LX + alpha * Lrho[: ,np.newaxis] * X
return fval, grad
~~~
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.
Then we call a solver to solve the nonlinear eigenvalue problem,
```python
# Execute the solver
X, out_dict = SLPG_smooth(obj_fun, M)
```
## Solvers
The `PySTOP` solver classes provide the solvers for optimization. Once we specify the Stiefel manifold and define the objective function, the `PySTOP` solver can be executed by
```
X, out_dict = name_of_solver(obj_fun, M)
```
Here `X` is the final output of the problem, and `out_dict` is a dictionary that contains the log information.
| Name | Comment | Call |
| ------------------------------------------------------------ | ------------------------------------------------------------ | ------------------ |
| [SLPG_smooth](https://stmopt.gitee.io/pystop/SLPG_smooth.html) | Penalty-free first-order method for smooth problems | `SLPG_smooth(...)` |
| [SLPG](https://stmopt.gitee.io/pystop/SLPG.html) | Penalty-free first-order method for nonsmooth problems | `SLPG(...)` |
| [SLPG_l21](https://stmopt.gitee.io/pystop/SLPG_l21.html) | Penalty-free first-order method for $\ell_{2,1}$-norm regularized problems | `SLPG_l21(...)` |
It is worth mentioning that `SLPG` solver supports customized nonsmooth regularization terms, interested reader could refer to the [description page](https://stmopt.gitee.io/pystop/SLPG.html) for details.
## Defining the objective function
Usually, computing the function value and gradient simultaneously is much faster than compute them separately, even when cache techniques are involved. Therefore, to achieve a better numerical performance, the existing solvers in PySTOP package requires an integrated call for function value and gradient of the objective function, i.e.
```python
# Defin objective function
def obj_fun(X):
'''
fval refers to the function value of f
grad refers to the gradient of f
'''
return fval, grad
```
Currently, PySTOP package does not involve build-in autodiff packages. However, we provides several useful functions in `pystop.utility` to help run the solvers only with specified objective function $f(X)$.
If you already know how to use NumPy, then it could be easy to use `autograd` package to generate the function value and gradient simultaneously. Just import `autograd.numpy` and setup the objective function by the build-in function provided in `autograd.numpy` to perform the computation. Once the function value is specified as `obj()`, we could apply the `pystop.utility.fun_autodiff` to generate the function that returns fval and grad simultaneously.
However, it is worth mentioning that `autograd.numpy` package only supports a subset of the standard NumPy package. Besides, the Autograd package does not support SciPy package. To achieve better numerical performance, we suggest the users to specify the `obj_fun` function manually.
```python
import autograd.numpy as anp
import numpy as np
n = 1000
p = 30
Z = np.random.randn(n, p)
def obj_fun(X):
return anp.sum((X-Z) **2 )
from pystop.utility import fun_autodiff
obj_grad, obj_fungrad = fun_autodiff(obj_fun)
from pystop.manifold import Stiefel
from pystop.solver import SLPG_smooth
M = Stiefel(n,p)
X, out_dict = SLPG_smooth(obj_fungrad, M)
```
Raw data
{
"_id": null,
"home_page": "https://stmopt.gitee.io/",
"name": "pystop",
"maintainer": "",
"docs_url": null,
"requires_python": ">=3.6",
"maintainer_email": "",
"keywords": "optimization,manifold optimization,Stiefel manifold",
"author": "Nachuan Xiao, Lei Wang, Bin Gao, Xin Liu, and Ya-xiang Yuan",
"author_email": "stmopt@foxmail.com",
"download_url": "https://files.pythonhosted.org/packages/2b/d6/3a001bd344d04e6e7f9aedf0be078c5284a68e39a7efcf3e69e17e15ff38/pystop-0.2.2.tar.gz",
"platform": null,
"description": "# PySTOP\n\n## Introduction\n\nThe STOP toolbox is designed for **optimization problems on the Stiefel manifold**, which could be expressed as \n$$\n\\begin{aligned}\n\t\\min_{X \\in \\mathbb{R}^{n\\times p}} ~ &f(X)\\\\\n\t\\text{s. t.}~& X^\\top X = I_p,\n\\end{aligned}\n$$\nwhere $I_p$ refers to the $p$-th order identity matrix, $X$ is a matrix with $n$ rows and $p$ columns. The feasible set of this optimization problem \n$$\n\\mathcal{S}_{n,p} := \\left\\{X \\in \\mathbb{R}^{n\\times p}: X^\\top X = I_p \\right\\},\n$$\ncan be regarded as a Riemannian manifold in $\\mathbb{R}^{n\\times p}$, and we also call it as **Stiefel manifold**. \n\nThis document describes the python version of the STOP package (PySTOP). Currently, PySTOP only involves the SLPG algorithm, which could handle smooth, $\\ell_1$-norm regularized and $\\ell_{2,1}$-nom regularized optimization problems on the Stiefel manifold. The package is \n\n## Installation\n\nThe source code of PySTOP package can be found from [the website](https://stmopt.gitee.io/). Besides, it supports direct installation from `pip`:\n\n```shell\npip install pystop\n```\n\n\n\n\n\n\n\n## Example\n\n### Problem formulation\n\nIn this section, we consider the following nonlinear eigenvalue problem\n$$\n\\min_{X \\in \\mathcal{S}_{n, p}} ~ \\frac{1}{2}\\mathrm{tr}(X^\\top L X) + \\frac{\\alpha}{4} \\rho^\\top L^{\\dagger} \\rho,\n$$\nwhere $\\rho = \\mathrm{Diag}(XX^\\top)$, and $L^{\\dagger}$ denotes the pseudo-inverse of the positive definite matrix $L$, i.e. $L^{\\dagger}LL^{\\dagger} = L^{\\dagger}$, $LL^{\\dagger}L = L$. Here we uses $\\mathrm{Diag}(M)$ to denote the vector that is composed of diagonal entries of the square matrix $M$, while $\\mathrm{diag}(v)$ refers to a diagonal matrix with $v$ to be its diagonal entries. Then the cost function and its **Euclidean gradient** can be expressed as\n$$\n\\begin{aligned}\n\t f(X) ={}& \\frac{1}{2}\\mathrm{tr}(X^\\top L X) + \\frac{\\alpha}{4} \\rho^\\top L^{\\dagger} \\rho,\\\\\n\t\\nabla f(X) ={}& LX + \\alpha \\mathrm{diag}(L^{\\dagger}\\rho)X.\n\\end{aligned}\n$$\n\n\nIn this example, we choose $L$ as a tri-diagonal matrix generated by `L = gallery('tridiag',n,-1,2,-1)`. Noting that $L$ is full-rank, then we can conclude that $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. For additional theory, readers are recommended to refer the papers in the about page. \n\n```python\n# Import packages \nimport numpy as np\nimport scipy as sp\nfrom scipy.sparse import diags\nfrom scipy.sparse.linalg import spsolve\n\n# Import manifolds and solvers\nfrom pystop.manifold import Stiefel\nfrom pystop.solver import SLPG_smooth\n\n\n# Set parameters\nn = 1000\np = 10\nalpha = 1\nM = Stiefel(n,p)\n\n# Defining objective function\nL = diags(np.array([-1, 2, -1]), np.array([1, 0, -1]), shape = (n,n)).tocsc()\ndef obj_fun(X):\n LX = L@X\n rho = np.sum(X * X, 1)\n Lrho = spsolve(L, rho)\n fval = 0.5*np.sum(X* LX) + (alpha /4) * np.sum(rho * Lrho)\n grad = LX + alpha * Lrho[: ,np.newaxis] * X\n return fval, grad\n\n# Execute the solver\nX, out_dict = SLPG_smooth(obj_fun, M)\n\n```\n\n\n\nLet us review the code step by step. First, we specify the dimension of the problem and specify the Stiefel manifold. In `pySTOP` package, we need to specify the dimension of the Stiefel manifold before executing the solver. The Stiefel manifold should be specified as the *STOP manifold class*, for example, \n\n```python\n# Set parameters\nn = 1000\np = 10\nalpha = 1\n# Specify the Stiefel manifold\nM = Stiefel(n,p)\n```\n\nHere `pystop.manifold.stiefel` is a build-in function to specify the Stiefel manifold and hence provides essential tools for the algorithm. \n\nThen we generate the data (matrix $L$) for the optimization problem by the following code,\n\n~~~python\nL = diags(np.array([-1, 2, -1]), np.array([1, 0, -1]), shape = (n,n)).tocsc()\n~~~\n\nHere we utilize `SciPy.sparse` to create a sparse representation of $L$ . Therefore, in each step we could use the `scipy.sparse.linalg.spsolve` function to compute . \n\nThen we specify the cost function and its gradient in the following function\n\n~~~python\n# Defin objective function\ndef obj_fun(X):\n LX = L@X\n rho = np.sum(X * X, 1)\n Lrho = spsolve(L, rho)\n fval = 0.5*np.sum(X* LX) + (alpha /4) * np.sum(rho * Lrho)\n grad = LX + alpha * Lrho[: ,np.newaxis] * X\n return fval, grad\n~~~\n\nCurrently, 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. \n\nThen we call a solver to solve the nonlinear eigenvalue problem, \n\n```python\n# Execute the solver\nX, out_dict = SLPG_smooth(obj_fun, M)\n```\n\n\n\n\n\n## Solvers\n\nThe `PySTOP` solver classes provide the solvers for optimization. Once we specify the Stiefel manifold and define the objective function, the `PySTOP` solver can be executed by \n\n```\nX, out_dict = name_of_solver(obj_fun, M)\n```\n\nHere `X` is the final output of the problem, and `out_dict` is a dictionary that contains the log information. \n\n| Name | Comment | Call |\n| ------------------------------------------------------------ | ------------------------------------------------------------ | ------------------ |\n| [SLPG_smooth](https://stmopt.gitee.io/pystop/SLPG_smooth.html) | Penalty-free first-order method for smooth problems | `SLPG_smooth(...)` |\n| [SLPG](https://stmopt.gitee.io/pystop/SLPG.html) | Penalty-free first-order method for nonsmooth problems | `SLPG(...)` |\n| [SLPG_l21](https://stmopt.gitee.io/pystop/SLPG_l21.html) | Penalty-free first-order method for $\\ell_{2,1}$-norm regularized problems | `SLPG_l21(...)` |\n\nIt is worth mentioning that `SLPG` solver supports customized nonsmooth regularization terms, interested reader could refer to the [description page](https://stmopt.gitee.io/pystop/SLPG.html) for details. \n\n\n\n## Defining the objective function\n\nUsually, computing the function value and gradient simultaneously is much faster than compute them separately, even when cache techniques are involved. Therefore, to achieve a better numerical performance, the existing solvers in PySTOP package requires an integrated call for function value and gradient of the objective function, i.e. \n\n```python\n# Defin objective function\ndef obj_fun(X):\n '''\n \tfval refers to the function value of f\n \tgrad refers to the gradient of f\n '''\n return fval, grad\n```\n\n\n\nCurrently, PySTOP package does not involve build-in autodiff packages. However, we provides several useful functions in `pystop.utility` to help run the solvers only with specified objective function $f(X)$. \n\nIf you already know how to use NumPy, then it could be easy to use `autograd` package to generate the function value and gradient simultaneously. Just import `autograd.numpy` and setup the objective function by the build-in function provided in `autograd.numpy` to perform the computation. Once the function value is specified as `obj()`, we could apply the `pystop.utility.fun_autodiff` to generate the function that returns fval and grad simultaneously.\n\nHowever, it is worth mentioning that `autograd.numpy` package only supports a subset of the standard NumPy package. Besides, the Autograd package does not support SciPy package. To achieve better numerical performance, we suggest the users to specify the `obj_fun` function manually. \n\n```python\nimport autograd.numpy as anp\nimport numpy as np\n\nn = 1000\np = 30\nZ = np.random.randn(n, p)\n\ndef obj_fun(X):\n return anp.sum((X-Z) **2 )\n\n\nfrom pystop.utility import fun_autodiff\n\nobj_grad, obj_fungrad = fun_autodiff(obj_fun)\n\n\nfrom pystop.manifold import Stiefel\nfrom pystop.solver import SLPG_smooth\n\nM = Stiefel(n,p)\nX, out_dict = SLPG_smooth(obj_fungrad, M)\n```\n\n\n\n\n\n\n\n",
"bugtrack_url": null,
"license": "",
"summary": "A Toolbox for Stiefel Manifold Optimization",
"version": "0.2.2",
"split_keywords": [
"optimization",
"manifold optimization",
"stiefel manifold"
],
"urls": [
{
"comment_text": "",
"digests": {
"md5": "b85fa455c9d39c8f57174f52436f4f55",
"sha256": "cd9465a1673261fb78944d3041581cfda8baf4628bfc7eb8545e46dcc0e6a011"
},
"downloads": -1,
"filename": "pystop-0.2.2-py3-none-any.whl",
"has_sig": false,
"md5_digest": "b85fa455c9d39c8f57174f52436f4f55",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": ">=3.6",
"size": 28972,
"upload_time": "2022-12-18T05:22:28",
"upload_time_iso_8601": "2022-12-18T05:22:28.540483Z",
"url": "https://files.pythonhosted.org/packages/7b/d4/6cb3b5e40b4375edbd2cba0548d1817664d4fba0589aad02547f8ff23e08/pystop-0.2.2-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
},
{
"comment_text": "",
"digests": {
"md5": "af1c1560d1a6ca7a98f4607194e1d320",
"sha256": "4087cf5140f5339f37691bcc6707c98cef8ecdfbce1b79a824c5f7090d708470"
},
"downloads": -1,
"filename": "pystop-0.2.2.tar.gz",
"has_sig": false,
"md5_digest": "af1c1560d1a6ca7a98f4607194e1d320",
"packagetype": "sdist",
"python_version": "source",
"requires_python": ">=3.6",
"size": 23850,
"upload_time": "2022-12-18T05:22:30",
"upload_time_iso_8601": "2022-12-18T05:22:30.440642Z",
"url": "https://files.pythonhosted.org/packages/2b/d6/3a001bd344d04e6e7f9aedf0be078c5284a68e39a7efcf3e69e17e15ff38/pystop-0.2.2.tar.gz",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2022-12-18 05:22:30",
"github": false,
"gitlab": false,
"bitbucket": false,
"lcname": "pystop"
}
JSON