diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index b271fb7..79968bf 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -72,8 +72,13 @@ jobs: echo "$(brew --prefix gcc)/bin" >> "$GITHUB_PATH" "$(brew --prefix gcc)/bin/gfortran" --version + # Without python-version every job resolves to whichever + # interpreter uv picks by default, so the matrix axis above was + # silently testing one version on all of them. - name: Setup uv uses: astral-sh/setup-uv@v7 + with: + python-version: ${{ matrix.python-version }} - name: Install project run: uv run bin/build.py install diff --git a/README.md b/README.md index 98704b2..d3f8b15 100644 --- a/README.md +++ b/README.md @@ -1,42 +1,115 @@ -# +![smopt](https://raw.githubusercontent.com/eggzec/smopt/master/docs/assets/images/smopt-banner.png) -**** +# smopt -[![Tests](https://github.com/eggzec//actions/workflows/test.yml/badge.svg)](https://github.com/eggzec//actions/workflows/test.yml) -[![Documentation](https://github.com/eggzec//actions/workflows/docs.yml/badge.svg)](https://github.com/eggzec//actions/workflows/docs.yml) +**Stiefel manifold optimization, with all numerics in Fortran 77** + +[![Tests](https://github.com/eggzec/smopt/actions/workflows/test.yml/badge.svg)](https://github.com/eggzec/smopt/actions/workflows/test.yml) +[![Documentation](https://github.com/eggzec/smopt/actions/workflows/docs.yml/badge.svg)](https://github.com/eggzec/smopt/actions/workflows/docs.yml) [![Ruff](https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/astral-sh/ruff/main/assets/badge/v2.json)](https://github.com/astral-sh/ruff) -[![codecov](https://codecov.io/github/eggzec//graph/badge.svg)](https://codecov.io/github/eggzec/) -[![Quality Gate Status](https://sonarcloud.io/api/project_badges/measure?project=eggzec_&metric=alert_status)](https://sonarcloud.io/project/overview?id=eggzec_) +[![codecov](https://codecov.io/github/eggzec/smopt/graph/badge.svg)](https://codecov.io/github/eggzec/smopt) +[![Quality Gate Status](https://sonarcloud.io/api/project_badges/measure?project=eggzec_smopt&metric=alert_status)](https://sonarcloud.io/project/overview?id=eggzec_smopt) [![License](https://img.shields.io/badge/license-GPL%203.0-blue.svg)](./LICENSE) -[![PyPI Downloads](https://img.shields.io/pypi/dm/.svg?label=PyPI%20downloads)](https://pypi.org/project//) -[![Python versions](https://img.shields.io/pypi/pyversions/.svg)](https://pypi.org/project//) +[![PyPI Downloads](https://img.shields.io/pypi/dm/smopt.svg?label=PyPI%20downloads)](https://pypi.org/project/smopt/) +[![Python versions](https://img.shields.io/pypi/pyversions/smopt.svg)](https://pypi.org/project/smopt/) + +`smopt` minimizes a smooth function, optionally plus a nonsmooth +regularizer, over the matrices with orthonormal columns: + +$$ +\min_{X \in \mathbb{R}^{n\times p}} f(X) + r(X) +\quad \text{subject to} \quad X^\top X = I_p. +$$ -`` +The solvers are penalty-free first-order methods. Rather than retracting +along geodesics, they work in the ambient space and dissolve the +orthogonality constraint with a cheap feasibility restoring map, so an +iteration costs little more than a gradient evaluation and a couple of +small matrix products. + +Everything numerical — the manifold geometry, the proximal operators, +the Barzilai-Borwein step sizes and the solver loops themselves — is +written in Fortran 77 and reached through f2py. Python supplies the +objective through a callback and handles reporting. NumPy is the only +runtime dependency; the extension links against nothing but the Fortran +runtime. ## Quick example ```python -import +import numpy as np +from smopt import slpg_smooth, Stiefel + +n, p = 1000, 10 +M = Stiefel(n, p) +A = np.diag(np.arange(n, dtype=float)) + + +def obj_fun(X): + """Return the objective and its Euclidean gradient together.""" + AX = A @ X + return float(np.sum(X * AX)), 2.0 * AX + + +X, out = slpg_smooth(obj_fun, M) +print(out["fval"], out["fea"]) ``` +## Solvers + +| Name | Comment | Call | +| --- | --- | --- | +| `slpg_smooth` | penalty-free first-order method for smooth problems | `slpg_smooth(obj_fun, M)` | +| `slpg` | penalty-free first-order method for nonsmooth problems | `slpg(obj_fun, M, prox=...)` | +| `slpg_l21` | penalty-free first-order method for $\ell_{2,1}$ regularized problems | `slpg_l21(obj_fun, M, gamma=...)` | +| `pencf` | constraint dissolving penalty method | `pencf(xinit, obj_fun, M)` | + +Every solver returns `(X, out)`, where `out` carries the `fvals`, `kkts` +and `feas` histories together with the final `fval`, `kkt` and `fea`. + ## Installation ```bash -pip install +pip install smopt ``` -Requires Python 3.10+ and NumPy. No external runtime dependencies. See the -[full installation guide](https://eggzec.github.io//installation/) for +Requires Python 3.10+ and NumPy. See the +[full installation guide](https://eggzec.github.io/smopt/installation/) for uv, poetry, and source builds. +Building from source additionally needs a Fortran compiler; the wheels +carry the Fortran runtime, so installing one does not. + +## Layout + +``` +src/ + smblas.f dense kernels: matmul, Cholesky, Jacobi eigensolver, Gram-Schmidt + smman.f Stiefel geometry: C, JA, JC, the A map, the polar retraction + smprox.f proximal operators and the l_{2,1} multiplier + smslpg.f the SLPG solver drivers and the Arrow-Hurwicz inner iteration + smpencf.f the pencf driver + _smopt.pyf f2py signatures binding the above to Python + smopt/ the thin Python layer: argument marshalling and reporting +tests/ + reference.py a NumPy transcription of the algorithm, used as a test oracle +``` + ## Documentation -- [Theory](https://eggzec.github.io//theory/) — mathematical background, hierarchical basis, algorithms -- [Quickstart](https://eggzec.github.io//quickstart/) — runnable examples -- [API Reference](https://eggzec.github.io//api/) — class and function signature and arguments -- [References](https://eggzec.github.io//references/) — literature citations +- [Theory](https://eggzec.github.io/smopt/theory/) — the manifold, the constraint dissolving map, the algorithms +- [Quickstart](https://eggzec.github.io/smopt/quickstart/) — runnable examples +- [API Reference](https://eggzec.github.io/smopt/api/) — class and function signatures and arguments +- [References](https://eggzec.github.io/smopt/references/) — literature citations + +## Acknowledgement + +The algorithm ported here originates in the STOP toolbox by Nachuan +Xiao, Lei Wang, Bin Gao, Xin Liu and Ya-xiang Yuan +(). `smopt` re-implements its numerics in +Fortran 77 behind the same solver interface. ## License diff --git a/bin/build.py b/bin/build.py index 83be46e..7437a14 100755 --- a/bin/build.py +++ b/bin/build.py @@ -78,7 +78,7 @@ def wheel(): def clean(): logger.debug("Starting cleanup ...") - run_command("uv pip uninstall ") + run_command("uv pip uninstall smopt") for entry in Path("").iterdir(): if entry.name in ["dist", "build", "lib", ".pytest_cache", ".ruff_cache"]: @@ -103,7 +103,7 @@ def clean(): def main(): - parser = argparse.ArgumentParser(description=" Build Script") + parser = argparse.ArgumentParser(description="smopt Build Script") parser.add_argument( "mode", help="""Build mode: diff --git a/bin/run_f2py.py b/bin/run_f2py.py new file mode 100644 index 0000000..81d9582 --- /dev/null +++ b/bin/run_f2py.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python3 +"""Generate the f2py wrappers for a signature file. + +f2py re-splits its own command line on whitespace, so a build directory +containing a space -- ``C:\\Users\\First Last\\...``, which is the norm on +Windows -- is torn into fragments and the run fails with a confusing +mixture of "Skipping file" and "Access is denied" errors. + +The workaround is to make sure no argument f2py sees ever contains a +separator or a space: the signature file is copied into the build +directory and f2py is invoked from there on the bare file name. + +Usage: + run_f2py.py [extra f2py arguments...] +""" + +import shutil +import subprocess +import sys +from pathlib import Path + + +def main() -> int: + """Copy the signature file into the build directory and run f2py. + + Returns: + The exit status of the f2py invocation. + """ + if len(sys.argv) < 3: + print(__doc__, file=sys.stderr) + return 2 + + signature = Path(sys.argv[1]).resolve() + build_dir = Path(sys.argv[2]).resolve() + extra_args = sys.argv[3:] + + build_dir.mkdir(parents=True, exist_ok=True) + staged = build_dir / signature.name + if staged != signature: + shutil.copyfile(signature, staged) + + return subprocess.call( # noqa: S603 + [ + sys.executable, + "-m", + "numpy.f2py", + signature.name, + "--build-dir", + ".", + *extra_args, + ], + cwd=str(build_dir), + ) + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/docs/api.md b/docs/api.md index b0d5c88..af878d8 100644 --- a/docs/api.md +++ b/docs/api.md @@ -1 +1,161 @@ # API Reference + +Everything below is re-exported at the top level, so +`from smopt import slpg_smooth` and +`from smopt.solver import slpg_smooth` are equivalent. + +## Manifold + +### `smopt.Stiefel` + +```python +Stiefel(n: int, p: int) +``` + +The manifold $\{X \in \mathbb{R}^{n\times p} : X^\top X = I_p\}$. Carries +the dimensions and exposes the geometric maps; each one is evaluated by +the Fortran 77 core. Raises `ValueError` if the dimensions are not +positive or if `p > n`. + +| Attribute | Meaning | +| --- | --- | +| `dim` | `n * p`, the dimension of the ambient space | + +| Method | Returns | +| --- | --- | +| `phi(m)` | $(M + M^\top)/2$ for a `(p, p)` matrix | +| `c(x)` | $X^\top X - I_p$ | +| `feas_eval(x)` | $\|X^\top X - I_p\|_F$, as a `float` | +| `ja(x, g)` | $G - X\,\Phi(X^\top G)$ | +| `jc(x, lam)` | $X\,\Phi(\Lambda)$ | +| `jc_transpose(x, d)` | $\Phi(X^\top D)$ | +| `a(x)` | the feasibility restoring map | +| `post_process(x)` | the orthogonal polar factor $UV^\top$ | +| `init_point(xinit=None)` | a feasible starting point | + +`init_point` draws a standard normal matrix when `xinit` is omitted, and +orthonormalizes whatever it ends up with unless it is already feasible. + +## Solvers + +Every solver returns `(X, out)`, where `X` is the solution and `out` is a +dictionary: + +| Key | Meaning | +| --- | --- | +| `fvals` | objective value per iteration | +| `kkts` | stationarity measure per iteration | +| `feas` | feasibility measure per iteration | +| `fval`, `kkt`, `fea` | the final values | +| `beta` | the penalty used, `pencf` only | + +The objective is a single callable returning the value and the Euclidean +gradient together: + +```python +def obj_fun(X): # X has shape (n, p) + return fval, grad # grad has shape (n, p) +``` + +`verbosity` is `0` for silence, `1` for the convergence and +post-processing lines, and `2` to also print periodically. `maxit` must +be at least `1`. A `xinit` you supply is used exactly as given; only the +default one is drawn and orthonormalized. + +### `smopt.slpg_smooth` + +```python +slpg_smooth( + obj_fun, + manifold, + xinit=None, + maxit=100, + gtol=1e-5, + post_process=True, + verbosity=2, +) +``` + +For a smooth objective. Prints every 20th iteration at `verbosity=2`. + +### `smopt.slpg` + +```python +slpg( + obj_fun, + manifold, + xinit=None, + maxit=100, + prox=None, + gtol=1e-5, + post_process=True, + verbosity=2, +) +``` + +For `f(X) + r(X)` with `r` reached through `prox(X, eta)`, which must +minimize $\|Y - X\|_F^2/(2\eta) + r(Y)$. `prox` defaults to the identity, +recovering the smooth case. Prints every 50th iteration at `verbosity=2`. + +### `smopt.slpg_l21` + +```python +slpg_l21( + obj_fun, + manifold, + xinit=None, + maxit=100, + gamma=0, + gtol=1e-5, + post_process=True, + verbosity=2, +) +``` + +For $f(X) + \gamma\|X\|_{2,1}$, which induces row sparsity. Prints every +50th iteration at `verbosity=2`. + +### `smopt.pencf` + +```python +pencf( + xinit, + obj_fun, + manifold, + beta=None, + maxit=100, + gtol=1e-5, + post_process=True, + verbosity=2, +) +``` + +A constraint dissolving penalty method. Note that the starting point +comes **first**. `beta` defaults to $0.1\|\nabla f(X_0)\|_F$; the value +actually used is reported as `out["beta"]`. Prints every 20th iteration +at `verbosity=2`. + +## Proximal operators + +### `smopt.prox_l1` + +```python +prox_l1(x, eta, gamma=0) +``` + +Proximal operator of $\gamma\|X\|_1$: entrywise soft thresholding. + +### `smopt.prox_l21` + +```python +prox_l21(x, eta, gamma=0) +``` + +Proximal operator of $\gamma\|X\|_{2,1}$: shrinks whole rows towards the +origin. + +Both are shaped so they can be handed straight to `slpg`: + +```python +X, out = slpg(obj_fun, M, prox=lambda X, eta: prox_l1(X, eta, gamma=0.05)) +``` diff --git a/docs/assets/images/smopt-banner.png b/docs/assets/images/smopt-banner.png new file mode 100644 index 0000000..1afe752 Binary files /dev/null and b/docs/assets/images/smopt-banner.png differ diff --git a/docs/assets/images/smopt-banner.svg b/docs/assets/images/smopt-banner.svg new file mode 100644 index 0000000..42272d2 --- /dev/null +++ b/docs/assets/images/smopt-banner.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/docs/assets/images/smopt-icon.svg b/docs/assets/images/smopt-icon.svg new file mode 100644 index 0000000..8170c25 --- /dev/null +++ b/docs/assets/images/smopt-icon.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/docs/assets/javascripts/katex.js b/docs/assets/javascripts/katex.js new file mode 100644 index 0000000..f973c74 --- /dev/null +++ b/docs/assets/javascripts/katex.js @@ -0,0 +1,32 @@ +/* Render the maths that pymdownx.arithmatex leaves in the page. + * + * mkdocs.yml sets arithmatex to generic mode, so formulas arrive as plain + * delimited text and KaTeX has to be run over the body. Material loads + * pages without a full reload, so the work is redone on every navigation + * via the document$ observable; the plain DOMContentLoaded path is the + * fallback for when that observable is not present. */ + +function renderMath(root) { + if (typeof renderMathInElement !== "function") return; + renderMathInElement(root, { + delimiters: [ + { left: "$$", right: "$$", display: true }, + { left: "$", right: "$", display: false }, + { left: "\\(", right: "\\)", display: false }, + { left: "\\[", right: "\\]", display: true }, + ], + // A failed formula should show as source, not blow up the page. + throwOnError: false, + ignoredTags: ["script", "noscript", "style", "textarea", "pre", "code"], + }); +} + +if (typeof document$ !== "undefined") { + document$.subscribe(function () { + renderMath(document.body); + }); +} else { + document.addEventListener("DOMContentLoaded", function () { + renderMath(document.body); + }); +} diff --git a/docs/assets/stylesheets/extra.css b/docs/assets/stylesheets/extra.css new file mode 100644 index 0000000..3d67d3d --- /dev/null +++ b/docs/assets/stylesheets/extra.css @@ -0,0 +1,91 @@ +/* smopt brand styling. + * + * Palette: Space berries, the same four colours the logo is built from. + * The header takes the deep berry, links and hover states take the hot + * pink, and the orange is kept for the few places that need to stand + * apart from both. */ + +:root { + --smopt-pink: #fd3db5; + --smopt-pink-light: #ffb8dc; + --smopt-orange: #fb6a2c; + --smopt-berry: #8c1946; +} + +/* Light scheme -------------------------------------------------------- */ + +[data-md-color-scheme="default"] { + --md-primary-fg-color: var(--smopt-berry); + --md-primary-fg-color--light: #a8265a; + --md-primary-fg-color--dark: #6d1236; + --md-accent-fg-color: var(--smopt-pink); + --md-typeset-a-color: #b31a6e; +} + +/* Dark scheme --------------------------------------------------------- */ + +[data-md-color-scheme="slate"] { + --md-primary-fg-color: var(--smopt-berry); + --md-primary-fg-color--light: #a8265a; + --md-primary-fg-color--dark: #5c0f2e; + --md-accent-fg-color: var(--smopt-pink); + --md-typeset-a-color: var(--smopt-pink); +} + +/* The logo is a dark tile, so give it a little breathing room rather + than letting it butt against the header text. */ +.md-header__button.md-logo img, +.md-header__button.md-logo svg { + border-radius: 4px; + height: 1.6rem; + width: 1.6rem; +} + +/* Tables carry the API and solver summaries, so make them legible + rather than decorative. */ +.md-typeset table:not([class]) th { + background-color: var(--smopt-berry); + color: #ffffff; + font-weight: 600; +} + +.md-typeset table:not([class]) tr:hover { + background-color: rgba(253, 61, 181, 0.06); +} + +/* Inline code reads as a term, not as a warning. */ +.md-typeset code { + border-radius: 3px; +} + +/* Admonitions in the brand's orange, which is otherwise unused and so + still registers as a signal. */ +.md-typeset .admonition.note, +.md-typeset details.note { + border-color: var(--smopt-orange); +} + +.md-typeset .note > .admonition-title, +.md-typeset .note > summary { + background-color: rgba(251, 106, 44, 0.1); +} + +.md-typeset .note > .admonition-title::before, +.md-typeset .note > summary::before { + background-color: var(--smopt-orange); +} + +/* KaTeX blocks are wide; let them scroll instead of stretching the page. */ +.md-typeset .arithmatex { + overflow-x: auto; + overflow-y: hidden; + padding: 0.2rem 0; +} + +/* The banner leads the index page. */ +.smopt-banner { + border-radius: 6px; + display: block; + margin: 0 auto 1.5rem; + max-width: 100%; +} diff --git a/docs/index.md b/docs/index.md index 9477141..b31d154 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,15 +1,60 @@ -# +# smopt -**** +![smopt](https://raw.githubusercontent.com/eggzec/smopt/master/docs/assets/images/smopt-banner.png) - +**Stiefel manifold optimization, with all numerics in Fortran 77** + +`smopt` minimizes a smooth function, optionally plus a nonsmooth +regularizer, over the set of matrices with orthonormal columns: + +$$ +\min_{X \in \mathbb{R}^{n\times p}} f(X) + r(X) +\quad \text{subject to} \quad X^\top X = I_p. +$$ ## Overview +The solvers are **penalty-free first-order methods**. Rather than +retracting along geodesics, they work in the ambient space and dissolve +the orthogonality constraint with a cheap feasibility restoring map, so +an iteration costs little more than a gradient evaluation and a couple +of small matrix products. + +Everything numerical — the manifold geometry, the proximal operators, +the Barzilai-Borwein step sizes and the solver loops themselves — is +implemented in Fortran 77 and reached through f2py. Python supplies the +objective through a callback and handles reporting; NumPy is the only +runtime dependency. + +```python +import numpy as np +from smopt import slpg_smooth, Stiefel + +M = Stiefel(1000, 10) +A = np.diag(np.arange(1000, dtype=float)) + + +def obj_fun(X): + AX = A @ X + return float(np.sum(X * AX)), 2.0 * AX + + +X, out = slpg_smooth(obj_fun, M) +``` + +## Solvers + +| Name | Use when | +| --- | --- | +| `slpg_smooth` | the objective is smooth | +| `slpg` | there is a nonsmooth term with a known proximal operator | +| `slpg_l21` | the nonsmooth term is $\gamma\|X\|_{2,1}$ | +| `pencf` | a constraint dissolving penalty method is preferred | + ## Documentation -- [Theory](theory.md) - mathematical background, hierarchical basis, algorithms +- [Theory](theory.md) - the manifold, the constraint dissolving map, the algorithms - [Installation](installation.md) - installation guide - [Quickstart](quickstart.md) - runnable examples -- [API Reference](api.md) - class and function signature and arguments +- [API Reference](api.md) - class and function signatures and arguments - [References](references.md) - literature citations diff --git a/docs/installation.md b/docs/installation.md index 39538b7..6bf258e 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -1,10 +1,10 @@ # Installation -`` can be installed from PyPI or directly from source via GitHub. +`smopt` can be installed from PyPI or directly from source via GitHub. --- -## [PyPI](https://pypi.org/project/) +## [PyPI](https://pypi.org/project/smopt) For using the PyPI package in your project, add it to your configuration file: @@ -12,7 +12,7 @@ For using the PyPI package in your project, add it to your configuration file: ```toml [project.dependencies] - = "*" # (1)! + smopt = "*" # (1)! ``` 1. Specifying a version is recommended @@ -20,7 +20,7 @@ For using the PyPI package in your project, add it to your configuration file: === "requirements.txt" ``` - >=0.1.0 + smopt>=0.1.0 ``` ### pip @@ -28,7 +28,7 @@ For using the PyPI package in your project, add it to your configuration file: === "Installation for user" ```bash - pip install --upgrade --user # (1)! + pip install --upgrade --user smopt # (1)! ``` 1. You may need to use `pip3` instead of `pip` depending on your Python installation. @@ -38,7 +38,7 @@ For using the PyPI package in your project, add it to your configuration file: ```bash python -m venv .venv source .venv/bin/activate - pip install --require-virtualenv --upgrade # (1)! + pip install --require-virtualenv --upgrade smopt # (1)! ``` 1. You may need to use `pip3` instead of `pip` depending on your Python installation. @@ -52,7 +52,7 @@ For using the PyPI package in your project, add it to your configuration file: === "Adding to uv project" ```bash - uv add + uv add smopt uv sync ``` @@ -60,41 +60,41 @@ For using the PyPI package in your project, add it to your configuration file: ```bash uv venv - uv pip install + uv pip install smopt ``` ### pipenv ```bash -pipenv install +pipenv install smopt ``` ### poetry ```bash -poetry add +poetry add smopt ``` ### pdm ```bash -pdm add +pdm add smopt ``` ### hatch ```bash -hatch add +hatch add smopt ``` --- -## [GitHub](https://github.com/eggzec/) +## [GitHub](https://github.com/eggzec/smopt) Install the latest development version directly from the repository: ```bash -pip install --upgrade "git+https://github.com/eggzec/.git#egg=" +pip install --upgrade "git+https://github.com/eggzec/smopt.git#egg=smopt" ``` ### Building locally @@ -102,8 +102,8 @@ pip install --upgrade "git+https://github.com/eggzec/.git#egg=.git -cd +git clone https://github.com/eggzec/smopt.git +cd smopt pip install -e . ``` diff --git a/docs/quickstart.md b/docs/quickstart.md index acb9843..f23f162 100644 --- a/docs/quickstart.md +++ b/docs/quickstart.md @@ -1 +1,162 @@ # Quickstart + +## A nonlinear eigenvalue problem + +The following problem is the standard smoke test for Stiefel solvers: + +$$ +\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}$ is the +pseudo-inverse of the positive definite $L$. The cost function and its +**Euclidean gradient** are + +$$ +\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} +$$ + +Taking $L$ tridiagonal, so that $L^{\dagger} = L^{-1}$: + +```python +import numpy as np +from scipy.sparse import diags +from scipy.sparse.linalg import spsolve + +from smopt import slpg_smooth, Stiefel + +n, p, alpha = 1000, 10, 1.0 +M = Stiefel(n, p) + +L = diags([-1, 2, -1], [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 + + +X, out = slpg_smooth(obj_fun, M) +``` + +!!! note + SciPy is only used to build this example's data; `smopt` itself + depends on NumPy alone. + +## Step by step + +### 1. Fix the manifold + +The dimensions are carried by a [`Stiefel`](api.md) instance, which also +exposes the geometric maps the solvers use: + +```python +from smopt import Stiefel + +M = Stiefel(1000, 10) +``` + +### 2. Define the objective + +A solver expects **one** callable returning the value and the Euclidean +gradient together: + +```python +def obj_fun(X): + ... + return fval, grad +``` + +Returning both at once is usually far cheaper than computing them +separately, even with caching, so `smopt` asks for them in a single +call. `X` arrives as an `(n, p)` array and `grad` must have the same +shape. + +### 3. Run a solver + +```python +X, out = slpg_smooth(obj_fun, M) +``` + +`X` is the solution and `out` is a dictionary of log information: + +| Key | Meaning | +| --- | --- | +| `fvals`, `kkts`, `feas` | per-iteration histories | +| `fval`, `kkt`, `fea` | the final objective, stationarity and feasibility | +| `beta` | the penalty `pencf` actually used | + +## Nonsmooth problems + +### A regularizer of your own + +Pass any `prox(X, eta)` that minimizes +$\|Y - X\|_F^2 / (2\eta) + r(Y)$. Here is $\ell_1$ regularization built +from the operator shipped with the package: + +```python +from smopt import prox_l1, slpg + +gamma = 0.05 +X, out = slpg(obj_fun, M, prox=lambda X, eta: prox_l1(X, eta, gamma=gamma)) +``` + +### Row sparsity + +For $r(X) = \gamma\|X\|_{2,1}$ use the dedicated driver, which knows the +prox and the constraint multiplier in closed form: + +```python +from smopt import slpg_l21 + +X, out = slpg_l21(obj_fun, M, gamma=1.0) +``` + +Whole rows of `X` are driven to zero, which selects variables: + +```python +import numpy as np + +live = np.linalg.norm(X, axis=1) > 1e-6 +print(f"{live.sum()} of {len(live)} rows survive") +``` + +## Choosing a solver + +| Solver | Use when | +| --- | --- | +| `slpg_smooth` | the objective is smooth | +| `slpg` | there is a nonsmooth term with a known prox | +| `slpg_l21` | the nonsmooth term is $\gamma\|X\|_{2,1}$ | +| `pencf` | a constraint dissolving penalty method is preferred | + +## Common options + +Every solver accepts: + +```python +X, out = slpg_smooth( + obj_fun, + M, + xinit=None, # starting point; random feasible point if omitted + maxit=100, # iteration budget + gtol=1e-5, # stationarity tolerance + post_process=True, # round the answer onto the manifold + verbosity=2, # 0 silent, 1 final lines, 2 periodic +) +``` + +`pencf` takes the starting point first and adds `beta`: + +```python +from smopt import pencf + +X, out = pencf(xinit, obj_fun, M, beta=None) +``` diff --git a/docs/references.md b/docs/references.md index b18be66..25cd276 100644 --- a/docs/references.md +++ b/docs/references.md @@ -1 +1,36 @@ # References + +The solvers implemented here follow the constraint dissolving and +penalty-free line of work on Stiefel manifold optimization. + +- Xiao, N., Liu, X., and Yuan, Y. (2022). *A class of smooth exact + penalty function methods for optimization problems with orthogonality + constraints.* Optimization Methods and Software, 37(4), 1205-1241. + +- Xiao, N., Liu, X., and Yuan, Y. (2022). *Exact penalty function for + $\ell_{2,1}$ norm minimization over the Stiefel manifold.* SIAM + Journal on Optimization, 31(4), 3097-3126. + +- Xiao, N., Liu, X., and Toh, K.-C. (2023). *Dissolving constraints for + Riemannian optimization.* Mathematics of Operations Research. + +- Barzilai, J. and Borwein, J. M. (1988). *Two-point step size gradient + methods.* IMA Journal of Numerical Analysis, 8(1), 141-148. + +- Edelman, A., Arias, T. A., and Smith, S. T. (1998). *The geometry of + algorithms with orthogonality constraints.* SIAM Journal on Matrix + Analysis and Applications, 20(2), 303-353. + +- Absil, P.-A., Mahony, R., and Sepulchre, R. (2008). *Optimization + Algorithms on Matrix Manifolds.* Princeton University Press. + +- Golub, G. H. and Van Loan, C. F. (2013). *Matrix Computations*, 4th + edition. Johns Hopkins University Press. Chapter 8 covers the cyclic + Jacobi eigenvalue algorithm used for the polar factor. + +## Provenance + +The algorithm ported here originates in the **STOP** toolbox by Nachuan +Xiao, Lei Wang, Bin Gao, Xin Liu and Ya-xiang Yuan, distributed at +. `smopt` re-implements its numerics in +Fortran 77 behind the same solver interface. diff --git a/docs/theory.md b/docs/theory.md index 006f9c5..5941b06 100644 --- a/docs/theory.md +++ b/docs/theory.md @@ -1 +1,134 @@ # Theory + +## The problem + +`smopt` solves + +$$ +\begin{aligned} + \min_{X \in \mathbb{R}^{n\times p}} ~ & f(X) + r(X)\ + \text{subject to}~& X^\top X = I_p, +\end{aligned} +$$ + +where $f$ is smooth and $r$ is a convex, possibly nonsmooth regularizer +reachable through its proximal operator. The feasible set + +$$ +\mathcal{S}_{n,p} := \left\{X \in \mathbb{R}^{n\times p}: X^\top X = I_p \right\} +$$ + +is the **Stiefel manifold**, a smooth embedded submanifold of +$\mathbb{R}^{n \times p}$ of dimension $np - p(p+1)/2$. + +## Constraint dissolving + +Classical Riemannian methods keep every iterate exactly on +$\mathcal{S}_{n,p}$, which costs a matrix decomposition per step. The +solvers here instead *dissolve* the constraint: they work in the ambient +space $\mathbb{R}^{n \times p}$ and rely on a cheap map that pulls a +drifting iterate back towards the manifold. + +That map is + +$$ +\mathcal{A}(X) = +\begin{cases} +\tfrac{3}{2} X - \tfrac{1}{2} X X^\top X, & \|X^\top X - I_p\|_F < \tfrac{1}{2},\[4pt] +X \left( \tfrac{1}{2}\left(X^\top X + I_p\right) \right)^{-1}, & \text{otherwise.} +\end{cases} +$$ + +The first branch is the second-order expansion of $X(X^\top X)^{-1/2}$ +about a feasible point and costs only matrix products. The second branch +is exact and needs a Cholesky factorization of order $p$, which is cheap +because $p \ll n$ in the problems of interest. Both branches fix the +manifold: $\mathcal{A}(X) = X$ whenever $X^\top X = I_p$. + +## Ingredients + +Writing $\Phi(M) = (M + M^\top)/2$ for the symmetrizing operator, the +solvers are built from + +| Map | Definition | Role | +| --- | --- | --- | +| $C(X)$ | $X^\top X - I_p$ | constraint violation | +| $\mathcal{J}_C(X)[\Lambda]$ | $X\,\Phi(\Lambda)$ | constraint Jacobian | +| $\mathcal{J}_C^\ast(X)[D]$ | $\Phi(X^\top D)$ | its adjoint | +| $\mathcal{J}_A(X)[G]$ | $G - X\,\Phi(X^\top G)$ | projected gradient | +| $\mathcal{A}(X)$ | above | feasibility restoration | + +Feasibility is measured by $\|C(X)\|_F$ throughout. + +## The solvers + +### SLPG + +The SLPG family takes a Barzilai-Borwein step along +$\mathcal{J}_A(X)[\nabla f(X)]$, applies the proximal operator of $r$, +and restores feasibility with $\mathcal{A}$. The step size uses one of +the two BB formulas + +$$ +\eta_k = \left| \frac{\langle S_k, Y_k\rangle}{\langle Y_k, Y_k \rangle} \right| +\qquad\text{or}\qquad +\eta_k = \left| \frac{\langle S_k, S_k\rangle}{\langle S_k, Y_k \rangle} \right|, +$$ + +with $S_k = X_k - X_{k-1}$ and $Y_k$ the corresponding change in the +search direction. The first few iterations use a conservative $c/L$ +instead, where $L$ is estimated from the gradient at the starting point. + +Three drivers are provided: + +- `slpg_smooth` for $r \equiv 0$. +- `slpg` for a general $r$, whose constraint multiplier $\Lambda$ is + tracked by an inner **Arrow-Hurwicz** iteration so that no penalty + parameter has to be tuned. +- `slpg_l21` for $r(X) = \gamma\|X\|_{2,1}$, where both the prox and the + multiplier + + $$ + \Lambda(X) = -\gamma\, X^\top \operatorname{diag}\!\left(\frac{1}{\|X_{i,:}\|_2}\right) X + $$ + + are available in closed form, so no inner iteration is needed. The + $\ell_{2,1}$ norm sums the Euclidean norms of the rows of $X$ and + therefore drives whole rows to zero, which is how sparse principal + component analysis and related models select variables. + +### pencf + +`pencf` adds an explicit penalty to the search direction, + +$$ +\mathcal{G}(X) = \mathcal{J}_A(X)[\nabla f(X)] + \beta\, \mathcal{J}_C(X)[C(X)], +$$ + +and restores feasibility only once $\|C(X)\|_F$ exceeds $10^{-1}$, +capping $\|X\|_F$ at $1.001\sqrt{p}$ to keep the iteration bounded. The +default $\beta$ is $0.1\|\nabla f(X_0)\|_F$. + +## Post-processing + +Because the iterates are only approximately feasible, every solver +optionally rounds the final point onto the manifold with the orthogonal +polar factor + +$$ +\mathcal{P}(X) = U V^\top = X \left(X^\top X\right)^{-1/2}, +\qquad X = U \Sigma V^\top, +$$ + +which is the nearest point of $\mathcal{S}_{n,p}$ in the Frobenius norm. +`smopt` computes it from a Jacobi eigendecomposition of the $p \times p$ +matrix $X^\top X$ rather than from a singular value decomposition of +$X$, which keeps the cost at $O(np^2 + p^3)$. + +## Implementation + +Every formula on this page is evaluated in Fortran 77, including the +solver loops themselves. Python supplies the objective through a +callback, and the linear algebra — matrix products, Cholesky, the Jacobi +eigensolver, modified Gram-Schmidt — is hand-written in the same +sources, so the extension links against nothing but the Fortran runtime. diff --git a/meson.build b/meson.build index 144a461..c51c048 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,11 @@ project( - '', + 'smopt', 'fortran', 'c', version: run_command(['python', 'bin/get_version.py'], check: true).stdout().strip() ) +# Passed down to src/ ; see the script for the f2py quoting bug it works +# around. +f2py_runner = files('bin/run_f2py.py') + subdir('src') diff --git a/mkdocs.yml b/mkdocs.yml index 0623e93..19ff5cd 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -1,11 +1,11 @@ -site_name: -site_description: -site_url: https://eggzec.github.io// +site_name: smopt +site_description: Stiefel manifold optimization with all numerics in Fortran 77 +site_url: https://eggzec.github.io/smopt/ docs_dir: docs site_dir: site -repo_name: eggzec/ -repo_url: https://github.com/eggzec/ +repo_name: eggzec/smopt +repo_url: https://github.com/eggzec/smopt # Copyright copyright: Copyright © 2026 eggzec @@ -51,8 +51,8 @@ theme: name: Switch to light mode # Logo and favicon - favicon: assets/images/.ico - logo: assets/images/.png + favicon: assets/images/smopt-icon.svg + logo: assets/images/smopt-icon.svg # Font configuration font: @@ -101,9 +101,9 @@ extra_css: # https://squidfunk.github.io/mkdocs-material/customization/#addit extra: social: - icon: fontawesome/brands/github - link: https://github.com/eggzec/ + link: https://github.com/eggzec/smopt - icon: fontawesome/brands/python - link: https://pypi.org/project// + link: https://pypi.org/project/smopt/ extra_javascript: - assets/javascripts/katex.js diff --git a/pyproject.toml b/pyproject.toml index 51a6442..a9a6534 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,9 +8,9 @@ requires = [ build-backend = "mesonpy" [project] -name = "" +name = "smopt" dynamic = ["version"] -description = "" +description = "Stiefel manifold optimization with all numerics in Fortran 77" authors = [ { name = "Saud Zahir", email = "m.saud.zahir@gmail.com" }, ] @@ -21,13 +21,16 @@ maintainers = [ readme = "README.md" license = "GPL-3.0" keywords = [ - "python" + "optimization", + "manifold optimization", + "Stiefel manifold", + "fortran", + "f2py" ] classifiers = [ "Development Status :: 5 - Production/Stable", "Intended Audience :: Science/Research", "Intended Audience :: Developers", - "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", "Programming Language :: Python", "Programming Language :: C", "Programming Language :: Fortran", @@ -50,12 +53,12 @@ dependencies = [ requires-python = ">=3.10" [project.urls] -homepage = "https://eggzec.github.io//" -documentation = "https://eggzec.github.io//" -source = "https://github.com/eggzec/" -changelog = "https://eggzec.github.io//changelog/" -releasenotes = "https://github.com/eggzec//releases/latest" -issues = "https://github.com/eggzec//issues" +homepage = "https://eggzec.github.io/smopt/" +documentation = "https://eggzec.github.io/smopt/" +source = "https://github.com/eggzec/smopt" +changelog = "https://eggzec.github.io/smopt/changelog/" +releasenotes = "https://github.com/eggzec/smopt/releases/latest" +issues = "https://github.com/eggzec/smopt/issues" [dependency-groups] dev = [ @@ -79,12 +82,12 @@ test = [ source = "vcs" [tool.hatch.build.targets.wheel] -include = [ "/" ] +include = [ "/smopt" ] [tool.ruff] -# CI workflow definitions and the helper scripts they run are not +# CI workflow definitions and the build helper scripts they run are not # part of the published package, so they are out of scope here. -extend-exclude = [".github"] +extend-exclude = [".github", "bin"] line-length = 80 indent-width = 4 preview = true @@ -125,7 +128,40 @@ ignore = [] [tool.ruff.lint.per-file-ignores] "**/__init__.py" = ["non-empty-init-module"] -"tests/*.py" = [] +"src/smopt/**/*.py" = [ + # The subpackages are one installable unit, so they reach each other + # with relative imports rather than by naming the distribution. + "relative-imports", + # The solver signatures are fixed by the algorithm this package + # ports; their arity is not ours to trim. + "too-many-arguments", + "too-many-positional-arguments", + # noqa is the portable spelling, and the one every reader knows. + "noqa-comments", +] +"tests/*.py" = [ + # A test asserts with assert. + "assert", + # Tolerances, sizes and iteration counts read better inline. + "magic-value-comparison", + # Annotating every fixture and parametrized case adds no clarity. + "ANN", + # The suite mirrors a public API that is deliberately CamelCase, and + # tests/reference.py is a faithful transcription of the original. + "invalid-function-name", + "invalid-argument-name", + "non-lowercase-variable-in-function", + "docstring-missing-returns", + "too-many-arguments", + "too-many-locals", + "too-many-positional-arguments", + # Parametrizing over a boolean flag is the point of those tests. + "boolean-type-hint-positional-argument", + "boolean-default-value-positional-argument", + "no-self-use", + "non-augmented-assignment", + "compare-to-empty-string", +] [tool.ruff.format] docstring-code-format = true diff --git a/src/_smopt.pyf b/src/_smopt.pyf new file mode 100644 index 0000000..0b0e757 --- /dev/null +++ b/src/_smopt.pyf @@ -0,0 +1,280 @@ +! -*- f90 -*- +! f2py signature file for the SMOPT Fortran 77 core. +! +! Every workspace array is declared intent(cache,hide) so that it is +! allocated by f2py and never appears in the Python signature. The +! objective, the proximal operator and the progress logger are Python +! callbacks; matrices cross the boundary flattened in column major +! order, which is how the Fortran side stores them anyway. + +python module __user__routines + interface + subroutine objfun(np,x,f,g) + integer intent(in) :: np + double precision dimension(np),intent(in) :: x + double precision intent(out) :: f + double precision dimension(np),intent(out),depend(np) :: g + end subroutine objfun + subroutine proxfn(np,x,eta,y) + integer intent(in) :: np + double precision dimension(np),intent(in) :: x + double precision intent(in) :: eta + double precision dimension(np),intent(out),depend(np) :: y + end subroutine proxfn + subroutine logfun(it,fval,kkt,fea,stage) + integer intent(in) :: it + double precision intent(in) :: fval + double precision intent(in) :: kkt + double precision intent(in) :: fea + integer intent(in) :: stage + end subroutine logfun + end interface +end python module __user__routines + +python module _smopt + interface + + subroutine smsymm(p,a) + integer intent(hide),depend(a) :: p = shape(a,0) + double precision dimension(p,p),intent(in,out,copy) :: a + end subroutine smsymm + + subroutine smcmap(n,p,x,c) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision dimension(p,p),intent(out),depend(p) :: c + end subroutine smcmap + + double precision function smfeas(n,p,x,wp) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp + end function smfeas + + subroutine smja(n,p,x,g,r,wp) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision dimension(n,p),intent(in),depend(n,p) :: g + double precision dimension(n,p),intent(out),depend(n,p) :: r + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp + end subroutine smja + + subroutine smjc(n,p,x,lam,r,wp) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision dimension(p,p),intent(in),depend(p) :: lam + double precision dimension(n,p),intent(out),depend(n,p) :: r + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp + end subroutine smjc + + subroutine smjct(n,p,x,d,r) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision dimension(n,p),intent(in),depend(n,p) :: d + double precision dimension(p,p),intent(out),depend(p) :: r + end subroutine smjct + + subroutine smamap(n,p,x,wp1,wp2,row) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p),intent(cache,hide),depend(p) :: row + end subroutine smamap + + subroutine smfix(n,p,x,wp1,wp2,row) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p),intent(cache,hide),depend(p) :: row + end subroutine smfix + + subroutine smpost(n,p,x,wp1,wp2,wp3,w,row) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp3 + double precision dimension(p),intent(cache,hide),depend(p) :: w + double precision dimension(p),intent(cache,hide),depend(p) :: row + end subroutine smpost + + subroutine sminit(n,p,x,wp) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp + end subroutine sminit + + subroutine smpl1(n,p,x,eta,gam,y) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision intent(in) :: eta + double precision intent(in) :: gam + double precision dimension(n,p),intent(out),depend(n,p) :: y + end subroutine smpl1 + + subroutine smpl21(n,p,x,eta,gam,eps,y) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision intent(in) :: eta + double precision intent(in) :: gam + double precision intent(in) :: eps + double precision dimension(n,p),intent(out),depend(n,p) :: y + end subroutine smpl21 + + subroutine smlm21(n,p,x,gam,lam) + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in) :: x + double precision intent(in) :: gam + double precision dimension(p,p),intent(out),depend(p) :: lam + end subroutine smlm21 + + subroutine smslps(n,p,x,maxit,gtol,ipost,objfun,logfun,nit,fvals,kkts,feasv,fval,kkt,fea,gf,gr,grp,s,y,xp,wp1,wp2,wp3,weig,row) + use __user__routines + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + integer intent(in) :: maxit + double precision intent(in) :: gtol + integer intent(in) :: ipost + external objfun + external logfun + integer intent(out) :: nit + double precision dimension(maxit),intent(out),depend(maxit) :: fvals + double precision dimension(maxit),intent(out),depend(maxit) :: kkts + double precision dimension(maxit),intent(out),depend(maxit) :: feasv + double precision intent(out) :: fval + double precision intent(out) :: kkt + double precision intent(out) :: fea + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gf + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gr + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: grp + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: s + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: y + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: xp + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp3 + double precision dimension(p),intent(cache,hide),depend(p) :: weig + double precision dimension(p),intent(cache,hide),depend(p) :: row + end subroutine smslps + + subroutine smslpg(n,p,x,maxit,gtol,ipost,objfun,proxfn,logfun,nit,fvals,kkts,feasv,fval,kkt,fea,gf,gr,grad,grdp,s,y,xp,z,xt,dx,lam,wp1,wp2,wp3,weig,row,steps) + use __user__routines + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + integer intent(in) :: maxit + double precision intent(in) :: gtol + integer intent(in) :: ipost + external objfun + external proxfn + external logfun + integer intent(out) :: nit + double precision dimension(maxit),intent(out),depend(maxit) :: fvals + double precision dimension(maxit),intent(out),depend(maxit) :: kkts + double precision dimension(maxit),intent(out),depend(maxit) :: feasv + double precision intent(out) :: fval + double precision intent(out) :: kkt + double precision intent(out) :: fea + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gf + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gr + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: grad + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: grdp + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: s + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: y + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: xp + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: z + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: xt + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: dx + double precision dimension(p,p),intent(cache,hide),depend(p) :: lam + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp3 + double precision dimension(p),intent(cache,hide),depend(p) :: weig + double precision dimension(p),intent(cache,hide),depend(p) :: row + double precision dimension(maxit),intent(cache,hide),depend(maxit) :: steps + end subroutine smslpg + + subroutine smsl21(n,p,x,maxit,gam,gtol,ipost,objfun,logfun,nit,fvals,kkts,feasv,fval,kkt,fea,gf,gr,grad,grdp,s,y,xp,xt,lam,wp1,wp2,wp3,weig,row) + use __user__routines + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + integer intent(in) :: maxit + double precision intent(in) :: gam + double precision intent(in) :: gtol + integer intent(in) :: ipost + external objfun + external logfun + integer intent(out) :: nit + double precision dimension(maxit),intent(out),depend(maxit) :: fvals + double precision dimension(maxit),intent(out),depend(maxit) :: kkts + double precision dimension(maxit),intent(out),depend(maxit) :: feasv + double precision intent(out) :: fval + double precision intent(out) :: kkt + double precision intent(out) :: fea + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gf + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gr + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: grad + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: grdp + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: s + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: y + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: xp + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: xt + double precision dimension(p,p),intent(cache,hide),depend(p) :: lam + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp3 + double precision dimension(p),intent(cache,hide),depend(p) :: weig + double precision dimension(p),intent(cache,hide),depend(p) :: row + end subroutine smsl21 + + subroutine smpcf(n,p,x,beta,maxit,gtol,ipost,objfun,logfun,nit,fvals,kkts,feasv,fval,kkt,fea,betout,gf,gr,gc,grp,s,y,xp,wp1,wp2,wp3,weig,row) + use __user__routines + integer intent(hide),depend(x) :: n = shape(x,0) + integer intent(hide),depend(x) :: p = shape(x,1) + double precision dimension(n,p),intent(in,out,copy) :: x + double precision intent(in) :: beta + integer intent(in) :: maxit + double precision intent(in) :: gtol + integer intent(in) :: ipost + external objfun + external logfun + integer intent(out) :: nit + double precision dimension(maxit),intent(out),depend(maxit) :: fvals + double precision dimension(maxit),intent(out),depend(maxit) :: kkts + double precision dimension(maxit),intent(out),depend(maxit) :: feasv + double precision intent(out) :: fval + double precision intent(out) :: kkt + double precision intent(out) :: fea + double precision intent(out) :: betout + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gf + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gr + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: gc + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: grp + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: s + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: y + double precision dimension(n,p),intent(cache,hide),depend(n,p) :: xp + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp1 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp2 + double precision dimension(p,p),intent(cache,hide),depend(p) :: wp3 + double precision dimension(p),intent(cache,hide),depend(p) :: weig + double precision dimension(p),intent(cache,hide),depend(p) :: row + end subroutine smpcf + + end interface +end python module _smopt diff --git a/src/meson.build b/src/meson.build index 66d0296..7804c28 100644 --- a/src/meson.build +++ b/src/meson.build @@ -1,12 +1,9 @@ - add_languages('fortran') py_mod = import('python') py = py_mod.find_installation(pure: false) py_dep = py.dependency() -f2py = find_program('f2py', required: true) - incdir_numpy = run_command(py, ['-c', 'import numpy; print(numpy.get_include())'], check: true @@ -22,39 +19,64 @@ fortranobject_c = run_command(py, check: true ).stdout().strip() -py_mod_name = '' +py_pkg_name = 'smopt' +py_mod_name = '_smopt' -fortran_sources = run_command(py, ['-c', - 'import glob; print(" ".join(glob.glob("*.f")))' -], check: true).stdout().strip().split() -message('Fortran sources: ', fortran_sources) - -python_sources = run_command(py, ['-c', - 'import glob; print(" ".join(glob.glob("*.py")))' -], check: true).stdout().strip().split() -message('Python sources: ', python_sources) +# Listed explicitly rather than globbed so that Meson can track them as +# real build inputs and reconfigure when one of them changes. +fortran_sources = [ + 'smblas.f', + 'smman.f', + 'smprox.f', + 'smslpg.f', + 'smpencf.f', +] pyf_sources = [ - '.pyf' + py_mod_name + '.pyf' ] message('Fortran sources: ', fortran_sources + pyf_sources) +# Driven through bin/run_f2py.py rather than f2py directly: f2py splits +# its own command line on whitespace, so passing a build directory that +# contains a space fails outright. f2py_target = custom_target( py_mod_name + '_f2py', input: pyf_sources, output: [py_mod_name + 'module.c', py_mod_name + '-f2pywrappers.f'], - command: [f2py, '@INPUT@', '--lower', '--build-dir', meson.current_build_dir()] + command: [py, f2py_runner, '@INPUT@', meson.current_build_dir(), '--lower'] ) +# Windows wheels must carry the Fortran runtime: a machine with Python +# but no MinGW has no libgfortran, libquadmath or libwinpthread to load. +# Note that having them on PATH is not enough either -- CPython 3.8 and +# later do not search PATH when resolving an extension's dependencies -- +# so they have to be linked in, not merely present. +# +# The -static-lib* flags cover the libraries they name. libwinpthread +# needs more than that: libgfortran.a refers to it, and the gfortran +# driver appends its own -lwinpthread from the spec file after every +# argument Meson passes, at which point the linker is back in dynamic +# mode and picks up the DLL. Even plain -static failed to prevent that +# on the MinGW-Builds UCRT toolchain used by the CI runners. Pulling the +# archive in with --whole-archive defines those symbols unconditionally +# and independently of ordering, so the trailing -lwinpthread finds +# nothing left to import. +# +# These are gfortran driver flags, and Meson would otherwise link this +# mixed C/Fortran target with gcc, hence the pinned link language below. +fc = meson.get_compiler('fortran') + if host_machine.system() == 'windows' - fortran_link_args = [ + fortran_link_args = fc.get_supported_link_arguments([ '-static-libgfortran', '-static-libgcc', '-static-libquadmath', - '-Wl,-Bstatic', + ]) + [ + '-Wl,--whole-archive,-Bstatic', '-lwinpthread', - '-Wl,-Bdynamic', + '-Wl,--no-whole-archive,-Bdynamic', ] else fortran_link_args = [] @@ -65,6 +87,41 @@ py.extension_module( fortran_sources + [f2py_target, fortranobject_c], c_args: ['-I' + incdir_numpy, '-I' + incdir_f2py], link_args: fortran_link_args, + link_language: 'fortran', dependencies: py_dep, + subdir: py_pkg_name, install: true ) + +py.install_sources( + [ + 'smopt/__init__.py', + 'smopt/_bridge.py', + ], + subdir: py_pkg_name +) + +py.install_sources( + [ + 'smopt/manifold/__init__.py', + 'smopt/manifold/stiefel.py', + ], + subdir: py_pkg_name / 'manifold' +) + +py.install_sources( + [ + 'smopt/solver/__init__.py', + 'smopt/solver/pencf.py', + 'smopt/solver/slpg.py', + ], + subdir: py_pkg_name / 'solver' +) + +py.install_sources( + [ + 'smopt/utility/__init__.py', + 'smopt/utility/utility.py', + ], + subdir: py_pkg_name / 'utility' +) diff --git a/src/smblas.f b/src/smblas.f new file mode 100644 index 0000000..98751c1 --- /dev/null +++ b/src/smblas.f @@ -0,0 +1,330 @@ +c----------------------------------------------------------------------- +c Dense linear algebra kernels used throughout SMOPT. +c +c Everything the solvers need is implemented here, so the extension +c module builds against nothing but a Fortran compiler and keeps the +c dependency footprint of the project template. +c +c All matrices are double precision and stored column major, which +c is the layout f2py hands over from NumPy. +c----------------------------------------------------------------------- + + DOUBLE PRECISION FUNCTION SMFRO(N, M, A) +c Frobenius norm of the N by M matrix A. + INTEGER N, M + DOUBLE PRECISION A(N,M) + INTEGER I, J + DOUBLE PRECISION T + + T = 0.0D0 + DO 20 J = 1, M + DO 10 I = 1, N + T = T + A(I,J)*A(I,J) + 10 CONTINUE + 20 CONTINUE + SMFRO = DSQRT(T) + RETURN + END + + + DOUBLE PRECISION FUNCTION SMDOT(N, M, A, B) +c Sum of the elementwise product of two N by M matrices. + INTEGER N, M + DOUBLE PRECISION A(N,M), B(N,M) + INTEGER I, J + DOUBLE PRECISION T + + T = 0.0D0 + DO 20 J = 1, M + DO 10 I = 1, N + T = T + A(I,J)*B(I,J) + 10 CONTINUE + 20 CONTINUE + SMDOT = T + RETURN + END + + + SUBROUTINE SMCOPY(N, M, A, B) +c B := A for N by M matrices. + INTEGER N, M + DOUBLE PRECISION A(N,M), B(N,M) + INTEGER I, J + + DO 20 J = 1, M + DO 10 I = 1, N + B(I,J) = A(I,J) + 10 CONTINUE + 20 CONTINUE + RETURN + END + + + SUBROUTINE SMSCAL(N, M, A, S) +c A := S*A for the N by M matrix A. + INTEGER N, M + DOUBLE PRECISION A(N,M), S + INTEGER I, J + + DO 20 J = 1, M + DO 10 I = 1, N + A(I,J) = S*A(I,J) + 10 CONTINUE + 20 CONTINUE + RETURN + END + + + SUBROUTINE SMMM(N, K, M, A, B, C) +c C := A*B with A(N,K), B(K,M) and C(N,M). + INTEGER N, K, M + DOUBLE PRECISION A(N,K), B(K,M), C(N,M) + INTEGER I, J, L + DOUBLE PRECISION T + + DO 40 J = 1, M + DO 10 I = 1, N + C(I,J) = 0.0D0 + 10 CONTINUE + DO 30 L = 1, K + T = B(L,J) + IF (T .NE. 0.0D0) THEN + DO 20 I = 1, N + C(I,J) = C(I,J) + A(I,L)*T + 20 CONTINUE + END IF + 30 CONTINUE + 40 CONTINUE + RETURN + END + + + SUBROUTINE SMMTM(N, K, M, A, B, C) +c C := A**T*B with A(N,K), B(N,M) and C(K,M). + INTEGER N, K, M + DOUBLE PRECISION A(N,K), B(N,M), C(K,M) + INTEGER I, J, L + DOUBLE PRECISION T + + DO 30 J = 1, M + DO 20 I = 1, K + T = 0.0D0 + DO 10 L = 1, N + T = T + A(L,I)*B(L,J) + 10 CONTINUE + C(I,J) = T + 20 CONTINUE + 30 CONTINUE + RETURN + END + + + SUBROUTINE SMRMM(N, P, X, W, ROW) +c X := X*W in place, with X(N,P) and W(P,P). ROW is a workspace +c vector of length P holding one transformed row at a time, which +c keeps the update free of any N by P scratch storage. + INTEGER N, P + DOUBLE PRECISION X(N,P), W(P,P), ROW(P) + INTEGER I, J, K + DOUBLE PRECISION T + + DO 40 I = 1, N + DO 20 J = 1, P + T = 0.0D0 + DO 10 K = 1, P + T = T + X(I,K)*W(K,J) + 10 CONTINUE + ROW(J) = T + 20 CONTINUE + DO 30 J = 1, P + X(I,J) = ROW(J) + 30 CONTINUE + 40 CONTINUE + RETURN + END + + + SUBROUTINE SMSYMM(P, A) +c A := (A + A**T)/2, the symmetrizing operator written Phi in the +c accompanying papers. + INTEGER P + DOUBLE PRECISION A(P,P) + INTEGER I, J + DOUBLE PRECISION T + + DO 20 J = 1, P + DO 10 I = J+1, P + T = 0.5D0*(A(I,J) + A(J,I)) + A(I,J) = T + A(J,I) = T + 10 CONTINUE + 20 CONTINUE + RETURN + END + + + SUBROUTINE SMCHOL(P, A, INFO) +c Cholesky factorization A = L*L**T of a symmetric positive +c definite matrix. L overwrites the lower triangle of A. INFO is +c zero on success and the index of the failing column otherwise. + INTEGER P, INFO + DOUBLE PRECISION A(P,P) + INTEGER I, J, K + DOUBLE PRECISION T + + INFO = 0 + DO 40 J = 1, P + T = A(J,J) + DO 10 K = 1, J-1 + T = T - A(J,K)*A(J,K) + 10 CONTINUE + IF (T .LE. 0.0D0) THEN + INFO = J + RETURN + END IF + A(J,J) = DSQRT(T) + DO 30 I = J+1, P + T = A(I,J) + DO 20 K = 1, J-1 + T = T - A(I,K)*A(J,K) + 20 CONTINUE + A(I,J) = T/A(J,J) + 30 CONTINUE + 40 CONTINUE + RETURN + END + + + SUBROUTINE SMSOLR(P, L, N, X) +c Solve M*z = x for every row x of X(N,P), where the symmetric +c positive definite M is supplied through its Cholesky factor L. +c X is overwritten by the solutions, that is X := X*M**(-1). + INTEGER P, N + DOUBLE PRECISION L(P,P), X(N,P) + INTEGER I, J, K + DOUBLE PRECISION T + + DO 50 I = 1, N + DO 20 J = 1, P + T = X(I,J) + DO 10 K = 1, J-1 + T = T - L(J,K)*X(I,K) + 10 CONTINUE + X(I,J) = T/L(J,J) + 20 CONTINUE + DO 40 J = P, 1, -1 + T = X(I,J) + DO 30 K = J+1, P + T = T - L(K,J)*X(I,K) + 30 CONTINUE + X(I,J) = T/L(J,J) + 40 CONTINUE + 50 CONTINUE + RETURN + END + + + SUBROUTINE SMJACO(P, A, W, V) +c Cyclic Jacobi eigenvalue decomposition of the symmetric matrix A, +c producing A = V*diag(W)*V**T. A is destroyed. The matrices met +c here are of order P, the column count of the iterate, so a few +c sweeps always suffice. + INTEGER P + DOUBLE PRECISION A(P,P), W(P), V(P,P) + INTEGER I, J, K, ISW + DOUBLE PRECISION OFF, THETA, T, C, S, H, U1, U2 + DOUBLE PRECISION TOL + PARAMETER (TOL = 1.0D-30) + + DO 20 J = 1, P + DO 10 I = 1, P + V(I,J) = 0.0D0 + 10 CONTINUE + V(J,J) = 1.0D0 + 20 CONTINUE + + DO 100 ISW = 1, 100 + OFF = 0.0D0 + DO 40 J = 2, P + DO 30 I = 1, J-1 + OFF = OFF + A(I,J)*A(I,J) + 30 CONTINUE + 40 CONTINUE + IF (OFF .LE. TOL) GO TO 110 + + DO 90 J = 2, P + DO 80 I = 1, J-1 + IF (A(I,J) .NE. 0.0D0) THEN + H = A(J,J) - A(I,I) + IF (DABS(H) + DABS(A(I,J)) .EQ. DABS(H)) THEN + T = A(I,J)/H + ELSE + THETA = 0.5D0*H/A(I,J) + T = 1.0D0/(DABS(THETA) + $ + DSQRT(1.0D0 + THETA*THETA)) + IF (THETA .LT. 0.0D0) T = -T + END IF + C = 1.0D0/DSQRT(1.0D0 + T*T) + S = T*C + DO 50 K = 1, P + U1 = A(I,K) + U2 = A(J,K) + A(I,K) = C*U1 - S*U2 + A(J,K) = S*U1 + C*U2 + 50 CONTINUE + DO 60 K = 1, P + U1 = A(K,I) + U2 = A(K,J) + A(K,I) = C*U1 - S*U2 + A(K,J) = S*U1 + C*U2 + 60 CONTINUE + DO 70 K = 1, P + U1 = V(K,I) + U2 = V(K,J) + V(K,I) = C*U1 - S*U2 + V(K,J) = S*U1 + C*U2 + 70 CONTINUE + END IF + 80 CONTINUE + 90 CONTINUE + 100 CONTINUE + + 110 CONTINUE + DO 120 I = 1, P + W(I) = A(I,I) + 120 CONTINUE + RETURN + END + + + SUBROUTINE SMMGS(N, P, X) +c Modified Gram-Schmidt orthonormalization. X is overwritten by an +c N by P matrix with orthonormal columns. + INTEGER N, P + DOUBLE PRECISION X(N,P) + INTEGER I, J, K + DOUBLE PRECISION T + + DO 60 J = 1, P + DO 30 K = 1, J-1 + T = 0.0D0 + DO 10 I = 1, N + T = T + X(I,K)*X(I,J) + 10 CONTINUE + DO 20 I = 1, N + X(I,J) = X(I,J) - T*X(I,K) + 20 CONTINUE + 30 CONTINUE + T = 0.0D0 + DO 40 I = 1, N + T = T + X(I,J)*X(I,J) + 40 CONTINUE + T = DSQRT(T) + IF (T .GT. 0.0D0) THEN + DO 50 I = 1, N + X(I,J) = X(I,J)/T + 50 CONTINUE + END IF + 60 CONTINUE + RETURN + END diff --git a/src/smman.f b/src/smman.f new file mode 100644 index 0000000..4d96f66 --- /dev/null +++ b/src/smman.f @@ -0,0 +1,281 @@ +c----------------------------------------------------------------------- +c Geometry of the Stiefel manifold +c +c S(n,p) = { X in R^(n x p) : X**T*X = I_p }. +c +c The routines below are the Fortran 77 counterparts of the manifold +c class exposed by SMOPT: the constraint map C, its Jacobian JC and +c adjoint, the projection JA onto the tangent-like space, the +c feasibility restoring map A, and the polar retraction used to +c round the final iterate back onto the manifold. +c----------------------------------------------------------------------- + + SUBROUTINE SMCMAP(N, P, X, C) +c C := X**T*X - I, the constraint violation of X. + INTEGER N, P + DOUBLE PRECISION X(N,P), C(P,P) + INTEGER I + + CALL SMMTM(N, P, P, X, X, C) + DO 10 I = 1, P + C(I,I) = C(I,I) - 1.0D0 + 10 CONTINUE + RETURN + END + + + DOUBLE PRECISION FUNCTION SMFEAS(N, P, X, WP) +c Feasibility measure ||X**T*X - I||_F. WP is P by P workspace. + INTEGER N, P + DOUBLE PRECISION X(N,P), WP(P,P) + DOUBLE PRECISION SMFRO + EXTERNAL SMFRO + + CALL SMCMAP(N, P, X, WP) + SMFEAS = SMFRO(P, P, WP) + RETURN + END + + + SUBROUTINE SMJA(N, P, X, G, R, WP) +c R := G - X*Phi(X**T*G), the projection that turns a Euclidean +c gradient G at X into the search direction used by the solvers. + INTEGER N, P + DOUBLE PRECISION X(N,P), G(N,P), R(N,P), WP(P,P) + INTEGER I, J + + CALL SMMTM(N, P, P, X, G, WP) + CALL SMSYMM(P, WP) + CALL SMMM(N, P, P, X, WP, R) + DO 20 J = 1, P + DO 10 I = 1, N + R(I,J) = G(I,J) - R(I,J) + 10 CONTINUE + 20 CONTINUE + RETURN + END + + + SUBROUTINE SMJC(N, P, X, LAM, R, WP) +c R := X*Phi(LAM), the Jacobian of the constraint map applied to a +c multiplier LAM. LAM is left untouched. + INTEGER N, P + DOUBLE PRECISION X(N,P), LAM(P,P), R(N,P), WP(P,P) + + CALL SMCOPY(P, P, LAM, WP) + CALL SMSYMM(P, WP) + CALL SMMM(N, P, P, X, WP, R) + RETURN + END + + + SUBROUTINE SMJCT(N, P, X, D, R) +c R := Phi(X**T*D), the adjoint of SMJC applied to a direction D. + INTEGER N, P + DOUBLE PRECISION X(N,P), D(N,P), R(P,P) + + CALL SMMTM(N, P, P, X, D, R) + CALL SMSYMM(P, R) + RETURN + END + + + SUBROUTINE SMAMAP(N, P, X, WP1, WP2, ROW) +c X := A(X), the feasibility restoring map. Close to the manifold +c the cheap second order expansion 1.5*X - X*(X**T*X)/2 is used, +c and the exact map X*((X**T*X + I)/2)**(-1) otherwise. + INTEGER N, P + DOUBLE PRECISION X(N,P), WP1(P,P), WP2(P,P), ROW(P) + INTEGER I, J, INFO + DOUBLE PRECISION FEAS + DOUBLE PRECISION SMFRO + EXTERNAL SMFRO + + CALL SMMTM(N, P, P, X, X, WP1) + DO 20 J = 1, P + DO 10 I = 1, P + WP2(I,J) = WP1(I,J) + 10 CONTINUE + WP2(J,J) = WP2(J,J) - 1.0D0 + 20 CONTINUE + FEAS = SMFRO(P, P, WP2) + + IF (FEAS .LT. 0.5D0) THEN + DO 40 J = 1, P + DO 30 I = 1, P + WP2(I,J) = -0.5D0*WP1(I,J) + 30 CONTINUE + WP2(J,J) = WP2(J,J) + 1.5D0 + 40 CONTINUE + CALL SMRMM(N, P, X, WP2, ROW) + ELSE + DO 60 J = 1, P + DO 50 I = 1, P + WP2(I,J) = 0.5D0*WP1(I,J) + 50 CONTINUE + WP2(J,J) = WP2(J,J) + 0.5D0 + 60 CONTINUE + CALL SMCHOL(P, WP2, INFO) + IF (INFO .EQ. 0) CALL SMSOLR(P, WP2, N, X) + END IF + RETURN + END + + + SUBROUTINE SMFIX(N, P, X, WP1, WP2, ROW) +c Feasibility restoration used by PENCF. Unlike SMAMAP the map is +c applied only once X has drifted appreciably off the manifold, and +c the result is capped in Frobenius norm to keep the penalized +c iteration bounded. + INTEGER N, P + DOUBLE PRECISION X(N,P), WP1(P,P), WP2(P,P), ROW(P) + INTEGER I, J, INFO + DOUBLE PRECISION FEAS, XN, CAP + DOUBLE PRECISION SMFRO + EXTERNAL SMFRO + + CALL SMMTM(N, P, P, X, X, WP1) + DO 20 J = 1, P + DO 10 I = 1, P + WP2(I,J) = WP1(I,J) + 10 CONTINUE + WP2(J,J) = WP2(J,J) - 1.0D0 + 20 CONTINUE + FEAS = SMFRO(P, P, WP2) + + IF (FEAS .GT. 1.0D-1) THEN + IF (FEAS .LT. 0.5D0) THEN + DO 40 J = 1, P + DO 30 I = 1, P + WP2(I,J) = -0.5D0*WP1(I,J) + 30 CONTINUE + WP2(J,J) = WP2(J,J) + 1.5D0 + 40 CONTINUE + CALL SMRMM(N, P, X, WP2, ROW) + ELSE + DO 60 J = 1, P + DO 50 I = 1, P + WP2(I,J) = 0.5D0*WP1(I,J) + 50 CONTINUE + WP2(J,J) = WP2(J,J) + 0.5D0 + 60 CONTINUE + CALL SMCHOL(P, WP2, INFO) + IF (INFO .EQ. 0) CALL SMSOLR(P, WP2, N, X) + END IF + END IF + + CAP = 1.001D0*DSQRT(DBLE(P)) + XN = SMFRO(N, P, X) + IF (XN .GT. CAP) CALL SMSCAL(N, P, X, CAP/XN) + RETURN + END + + + SUBROUTINE SMPOST(N, P, X, WP1, WP2, WP3, W, ROW) +c X := U*V**T where X = U*S*V**T is a thin singular value +c decomposition. When X has full column rank that orthogonal polar +c factor equals X*(X**T*X)**(-1/2), obtained here from a Jacobi +c eigendecomposition of the P by P matrix X**T*X rather than from a +c decomposition of X itself. +c +c A rank deficient X has no polar factor: the directions belonging +c to a vanishing singular value are unconstrained. Those columns +c are filled with an arbitrary orthonormal completion so that the +c result still lands on the manifold, which is what a singular value +c decomposition would hand back as well. Regularized solvers reach +c this case whenever they zero out enough rows. + INTEGER N, P + DOUBLE PRECISION X(N,P), WP1(P,P), WP2(P,P), WP3(P,P) + DOUBLE PRECISION W(P), ROW(P) + INTEGER I, J, K, M + DOUBLE PRECISION T, WMAX, TOL, THR + DOUBLE PRECISION EPS + PARAMETER (EPS = 2.220446049250313D-16) + + CALL SMMTM(N, P, P, X, X, WP1) + CALL SMSYMM(P, WP1) + CALL SMJACO(P, WP1, W, WP2) + +c X := X*V, whose K-th column is the K-th singular value times the +c corresponding left singular vector. + CALL SMRMM(N, P, X, WP2, ROW) + + WMAX = 0.0D0 + DO 10 K = 1, P + IF (W(K) .GT. WMAX) WMAX = W(K) + 10 CONTINUE + TOL = DBLE(P)*EPS*WMAX + +c Normalize the columns that carry a nonzero singular value and +c blank the rest, so the completion below sees only real vectors. + DO 40 K = 1, P + IF (W(K) .GT. TOL) THEN + T = 1.0D0/DSQRT(W(K)) + DO 20 I = 1, N + X(I,K) = T*X(I,K) + 20 CONTINUE + ELSE + DO 30 I = 1, N + X(I,K) = 0.0D0 + 30 CONTINUE + END IF + 40 CONTINUE + +c Complete the blanked columns. Projecting the N coordinate axes +c off the P-1 columns already in place leaves a total squared norm +c of N-P+1, so at least one axis clears the threshold below. + THR = 0.5D0*DSQRT(DBLE(N-P+1)/DBLE(N)) + DO 130 K = 1, P + IF (W(K) .LE. TOL) THEN + DO 120 J = 1, N + DO 50 I = 1, N + X(I,K) = 0.0D0 + 50 CONTINUE + X(J,K) = 1.0D0 + DO 80 M = 1, P + IF (M .NE. K) THEN + T = 0.0D0 + DO 60 I = 1, N + T = T + X(I,M)*X(I,K) + 60 CONTINUE + DO 70 I = 1, N + X(I,K) = X(I,K) - T*X(I,M) + 70 CONTINUE + END IF + 80 CONTINUE + T = 0.0D0 + DO 90 I = 1, N + T = T + X(I,K)*X(I,K) + 90 CONTINUE + T = DSQRT(T) + IF (T .GT. THR) THEN + DO 100 I = 1, N + X(I,K) = X(I,K)/T + 100 CONTINUE + GO TO 130 + END IF + 120 CONTINUE + END IF + 130 CONTINUE + +c X := U*V**T. + DO 150 J = 1, P + DO 140 I = 1, P + WP3(I,J) = WP2(J,I) + 140 CONTINUE + 150 CONTINUE + CALL SMRMM(N, P, X, WP3, ROW) + RETURN + END + + + SUBROUTINE SMINIT(N, P, X, WP) +c Orthonormalize X unless it already sits on the Stiefel manifold. + INTEGER N, P + DOUBLE PRECISION X(N,P), WP(P,P) + DOUBLE PRECISION SMFEAS + EXTERNAL SMFEAS + + IF (SMFEAS(N, P, X, WP) .GT. 1.0D-6) CALL SMMGS(N, P, X) + RETURN + END diff --git a/src/smopt/__init__.py b/src/smopt/__init__.py new file mode 100644 index 0000000..53fc8f3 --- /dev/null +++ b/src/smopt/__init__.py @@ -0,0 +1,55 @@ +r"""SMOPT, a toolbox for optimization over the Stiefel manifold. + +The problems addressed here read + +.. math:: + + \min_{X \in \mathbb{R}^{n \times p}} f(X) + r(X) + \quad \text{subject to} \quad X^\top X = I_p, + +whose feasible set is the Stiefel manifold :math:`\mathcal{S}_{n,p}`. +The smooth part ``f`` is supplied as a Python callable returning the +function value and its Euclidean gradient together; the optional +nonsmooth part ``r`` is reached through its proximal operator. + +All of the numerics, including the solver iterations themselves, are +implemented in Fortran 77 and reached through f2py. Python is +responsible only for marshalling arguments, calling back into the user +objective, and reporting progress. + +Examples: + >>> import numpy as np + >>> from smopt import Stiefel, slpg_smooth + >>> manifold = Stiefel(6, 2) + >>> a = np.diag([5.0, 4.0, 3.0, 2.0, 1.0, 0.0]) + >>> def obj(x): + ... return float(np.sum(x * (a @ x))), 2.0 * (a @ x) + >>> x, out = slpg_smooth( + ... obj, manifold, xinit=np.arange(12.0).reshape(6, 2), verbosity=0 + ... ) + >>> bool(manifold.feas_eval(x) < 1e-8) + True +""" + +from importlib.metadata import PackageNotFoundError, version + +from .manifold import Stiefel +from .solver import pencf, slpg, slpg_l21, slpg_smooth +from .utility import prox_l1, prox_l21 + + +try: + __version__ = version("smopt") +except PackageNotFoundError: # pragma: no cover - source checkout + __version__ = "0.0.0" + +__all__: list[str] = [ + "Stiefel", + "__version__", + "pencf", + "prox_l1", + "prox_l21", + "slpg", + "slpg_l21", + "slpg_smooth", +] diff --git a/src/smopt/_bridge.py b/src/smopt/_bridge.py new file mode 100644 index 0000000..b1f2c25 --- /dev/null +++ b/src/smopt/_bridge.py @@ -0,0 +1,175 @@ +"""Plumbing between the Python API and the Fortran 77 core. + +The Fortran drivers reach back into Python for three things: the +objective, the proximal operator of a nonsmooth regularizer, and +progress reporting. Matrices cross that boundary flattened in column +major order, which is how Fortran stores them, so the adapters here +reshape in both directions and keep every user facing array in the +familiar ``(n, p)`` form. +""" + +from collections.abc import Callable + +import numpy as np +from numpy.typing import NDArray + + +Matrix = NDArray[np.float64] +ObjFun = Callable[[Matrix], tuple[float, Matrix]] +Prox = Callable[[Matrix, float], Matrix] + +_LINE = "Iter:{} fval:{:.3e} kkts:{:.3e} feas:{:3e}" + +#: ``stage`` values passed to the logging callback by the Fortran side. +_STAGE_ITERATION = 0 +_STAGE_CONVERGED = 1 + +#: Verbosity at which the periodic progress lines are printed. +_VERBOSE_PERIODIC = 2 + + +def as_matrix(x: Matrix, n: int, p: int, name: str = "X") -> Matrix: + """Return ``x`` as a contiguous ``(n, p)`` array of doubles. + + Args: + x: Array to validate and convert. + n: Expected row count. + p: Expected column count. + name: Name used in the error message. + + Returns: + A float64 array of shape ``(n, p)``. + + Raises: + ValueError: If ``x`` does not have shape ``(n, p)``. + """ + out = np.asarray(x, dtype=np.float64) + if out.shape != (n, p): + msg = f"{name} must have shape {(n, p)}, got {out.shape}" + raise ValueError(msg) + return out + + +def check_maxit(maxit: int) -> int: + """Validate the iteration budget. + + Args: + maxit: Requested maximum number of iterations. + + Returns: + The validated budget. + + Raises: + ValueError: If ``maxit`` is not at least one. + """ + maxit = int(maxit) + if maxit < 1: + msg = f"maxit must be at least 1, got {maxit}" + raise ValueError(msg) + return maxit + + +def obj_callback(obj_fun: ObjFun, n: int, p: int) -> Callable[..., object]: + """Adapt a user objective for the Fortran ``OBJFUN`` callback. + + Args: + obj_fun: Callable mapping ``X`` to ``(fval, grad)``. + n: Row count of the iterate. + p: Column count of the iterate. + + Returns: + A callable taking the flattened iterate and returning the + function value together with the flattened gradient. + """ + + def objfun(x: Matrix) -> tuple[float, Matrix]: + fval, grad = obj_fun(x.reshape((n, p), order="F")) + flat = np.asarray(grad, dtype=np.float64).reshape(-1, order="F") + return float(fval), flat + + return objfun + + +def prox_callback(prox: Prox, n: int, p: int) -> Callable[..., object]: + """Adapt a user proximal operator for the Fortran ``PROXFN`` callback. + + Args: + prox: Callable mapping ``(X, eta)`` to the proximal point. + n: Row count of the iterate. + p: Column count of the iterate. + + Returns: + A callable taking the flattened point and the step size, and + returning the flattened proximal point. + """ + + def proxfn(x: Matrix, eta: float) -> Matrix: + y = prox(x.reshape((n, p), order="F"), float(eta)) + return np.asarray(y, dtype=np.float64).reshape(-1, order="F") + + return proxfn + + +def log_callback(verbosity: int, period: int) -> Callable[..., object]: + """Build the progress reporter handed to the Fortran drivers. + + The Fortran side reports every iteration and leaves the decision of + what to print here, so the printing policy stays in Python. + + Args: + verbosity: ``0`` silences output, ``1`` prints only the final + and post-processing lines, ``2`` also prints periodically. + period: Iteration stride between periodic lines. + + Returns: + A callable accepting ``(it, fval, kkt, fea, stage)``. + """ + + def logfun( + it: int, fval: float, kkt: float, fea: float, stage: int + ) -> None: + if stage == _STAGE_ITERATION: + if verbosity == _VERBOSE_PERIODIC and it % period == 0: + print(_LINE.format(it, fval, kkt, fea)) + elif stage == _STAGE_CONVERGED: + if verbosity >= 1: + print(_LINE.format(it, fval, kkt, fea)) + elif verbosity >= 1: + print("Post-processing") + print(_LINE.format(it, fval, kkt, fea)) + + return logfun + + +def output_dict( + nit: int, + fvals: Matrix, + kkts: Matrix, + feasv: Matrix, + fval: float, + kkt: float, + fea: float, +) -> dict[str, object]: + """Assemble the log dictionary returned alongside the solution. + + Args: + nit: Number of iterations actually performed. + fvals: Objective value history, padded to the iteration budget. + kkts: Stationarity history, padded to the iteration budget. + feasv: Feasibility history, padded to the iteration budget. + fval: Final objective value. + kkt: Final stationarity measure. + fea: Final feasibility measure. + + Returns: + A dictionary with the per-iteration histories truncated to the + iterations that ran, plus the final scalars. + """ + return { + "kkts": kkts[:nit].tolist(), + "fvals": fvals[:nit].tolist(), + "fea": fea, + "kkt": kkt, + "fval": fval, + "feas": feasv[:nit].tolist(), + } diff --git a/src/smopt/manifold/__init__.py b/src/smopt/manifold/__init__.py new file mode 100644 index 0000000..cf4580a --- /dev/null +++ b/src/smopt/manifold/__init__.py @@ -0,0 +1,6 @@ +"""Manifolds on which SMOPT can optimize.""" + +from .stiefel import Stiefel + + +__all__: list[str] = ["Stiefel"] diff --git a/src/smopt/manifold/stiefel.py b/src/smopt/manifold/stiefel.py new file mode 100644 index 0000000..c1b5976 --- /dev/null +++ b/src/smopt/manifold/stiefel.py @@ -0,0 +1,165 @@ +"""The Stiefel manifold and the maps the solvers need on it.""" + +import numpy as np + +from .. import _smopt +from .._bridge import Matrix, as_matrix + + +class Stiefel: + r"""The Stiefel manifold :math:`\{X \in R^{n \times p} : X^T X = I_p\}`. + + The object carries the dimensions and exposes the maps the solvers + need. Every one of them is evaluated by the Fortran 77 core, and each + is named after the symbol it carries in the theory documentation. + + Args: + n: Number of rows of the iterate. + p: Number of columns of the iterate. + + Raises: + ValueError: If the dimensions are not positive or ``p`` exceeds + ``n``. + + Examples: + >>> import numpy as np + >>> from smopt import Stiefel + >>> manifold = Stiefel(4, 2) + >>> x = manifold.init_point(np.eye(4, 2)) + >>> bool(manifold.feas_eval(x) < 1e-12) + True + """ + + def __init__(self, n: int, p: int) -> None: + n, p = int(n), int(p) + if n < 1 or p < 1: + msg = f"n and p must be positive, got n={n}, p={p}" + raise ValueError(msg) + if p > n: + msg = f"p must not exceed n, got n={n}, p={p}" + raise ValueError(msg) + self._n = n + self._p = p + self.dim = n * p + + def _mat(self, x: Matrix, name: str = "x") -> Matrix: + return as_matrix(x, self._n, self._p, name) + + def _sq(self, m: Matrix, name: str) -> Matrix: + return as_matrix(m, self._p, self._p, name) + + def phi(self, m: Matrix) -> Matrix: + """Symmetrize a square matrix. + + Args: + m: A ``(p, p)`` matrix. + + Returns: + ``(m + m.T) / 2``. + """ + return _smopt.smsymm(self._sq(m, "m")) + + def a(self, x: Matrix) -> Matrix: + """Pull a point back towards the manifold. + + Close to the manifold a second order expansion is used, and the + exact map ``x ((x^T x + I) / 2)^-1`` otherwise. + + Args: + x: An ``(n, p)`` matrix. + + Returns: + The restored point. + """ + return _smopt.smamap(self._mat(x)) + + def ja(self, x: Matrix, g: Matrix) -> Matrix: + """Project a Euclidean gradient onto the search direction. + + Args: + x: An ``(n, p)`` matrix. + g: The Euclidean gradient at ``x``. + + Returns: + ``g - x phi(x^T g)``. + """ + return _smopt.smja(self._mat(x), self._mat(g, "g")) + + def jc(self, x: Matrix, lam: Matrix) -> Matrix: + """Apply the constraint Jacobian to a multiplier. + + Args: + x: An ``(n, p)`` matrix. + lam: A ``(p, p)`` multiplier. + + Returns: + ``x phi(lam)``. + """ + return _smopt.smjc(self._mat(x), self._sq(lam, "lam")) + + def jc_transpose(self, x: Matrix, d: Matrix) -> Matrix: + """Apply the adjoint of the constraint Jacobian to a direction. + + Args: + x: An ``(n, p)`` matrix. + d: An ``(n, p)`` direction. + + Returns: + ``phi(x^T d)``. + """ + return _smopt.smjct(self._mat(x), self._mat(d, "d")) + + def c(self, x: Matrix) -> Matrix: + """Evaluate the constraint violation. + + Args: + x: An ``(n, p)`` matrix. + + Returns: + ``x^T x - I``. + """ + return _smopt.smcmap(self._mat(x)) + + def feas_eval(self, x: Matrix) -> float: + """Measure how far a point sits from the manifold. + + Args: + x: An ``(n, p)`` matrix. + + Returns: + The Frobenius norm of ``x^T x - I``. + """ + return float(_smopt.smfeas(self._mat(x))) + + def init_point(self, xinit: Matrix | None = None) -> Matrix: + """Produce a feasible starting point. + + Args: + xinit: Optional starting matrix. A standard normal matrix is + drawn when it is omitted. Either way the result is + orthonormalized unless it is already feasible. + + Returns: + An ``(n, p)`` matrix on the manifold. + """ + start = ( + np.random.randn(self._n, self._p) + if xinit is None + else self._mat(xinit, "xinit") + ) + return _smopt.sminit(start) + + def post_process(self, x: Matrix) -> Matrix: + """Round a point onto the manifold. + + Args: + x: An ``(n, p)`` matrix. + + Returns: + The orthogonal polar factor ``u v^T`` of ``x``, where + ``x = u s v^T`` is a thin singular value decomposition. + """ + return _smopt.smpost(self._mat(x)) + + +__all__: list[str] = ["Stiefel"] diff --git a/src/smopt/solver/__init__.py b/src/smopt/solver/__init__.py new file mode 100644 index 0000000..2df925e --- /dev/null +++ b/src/smopt/solver/__init__.py @@ -0,0 +1,7 @@ +"""Solvers for optimization over the Stiefel manifold.""" + +from .pencf import pencf +from .slpg import slpg, slpg_l21, slpg_smooth + + +__all__: list[str] = ["pencf", "slpg", "slpg_l21", "slpg_smooth"] diff --git a/src/smopt/solver/pencf.py b/src/smopt/solver/pencf.py new file mode 100644 index 0000000..d64f85a --- /dev/null +++ b/src/smopt/solver/pencf.py @@ -0,0 +1,100 @@ +"""PenCF, the constraint dissolving penalty solver.""" + +from typing import Any + +from .. import _smopt +from .._bridge import ( + Matrix, + ObjFun, + as_matrix, + check_maxit, + log_callback, + obj_callback, + output_dict, +) +from ..manifold import Stiefel + + +#: Iteration stride between the periodic progress lines. +_PERIOD = 20 + +#: Sentinel handed to the Fortran driver to request the default penalty. +_AUTO_BETA = -1.0 + + +def pencf( + xinit: Matrix, + obj_fun: ObjFun, + manifold: Stiefel, + beta: float | None = None, + maxit: int = 100, + gtol: float = 1e-5, + post_process: bool = True, # noqa: FBT001, FBT002 + verbosity: int = 2, + **kwargs: Any, # noqa: ANN401 +) -> tuple[Matrix, dict[str, object]]: + r"""Minimize a smooth objective with a constraint dissolving penalty. + + The search direction adds ``beta jc(x, c(x))`` to the projected + gradient, and feasibility is restored only once the iterate has + drifted appreciably off the manifold. + + Args: + xinit: Starting point. A random feasible point is drawn when it + is ``None``. + obj_fun: Callable mapping ``x`` to ``(fval, grad)``, where + ``grad`` is the Euclidean gradient. + manifold: The :class:`~smopt.manifold.Stiefel` instance fixing + the dimensions. + beta: Penalty weight. Defaults to ``0.1`` times the Frobenius + norm of the gradient at the starting point. + maxit: Maximum number of iterations. + gtol: Stationarity tolerance that stops the iteration. + post_process: Whether to round the final iterate onto the + manifold. + verbosity: ``0`` silences output, ``1`` prints the final lines, + ``2`` also prints every twentieth iteration. + **kwargs: Ignored, accepted so solvers stay interchangeable. + + Returns: + The solution and a dictionary holding the ``fvals``, ``kkts`` + and ``feas`` histories, the final ``fval``, ``kkt`` and ``fea`` + values, and the ``beta`` actually used. + + Examples: + >>> import numpy as np + >>> from smopt import Stiefel, pencf + >>> manifold = Stiefel(6, 2) + >>> a = np.diag([5.0, 4.0, 3.0, 2.0, 1.0, 0.0]) + >>> def obj(x): + ... return float(np.sum(x * (a @ x))), 2.0 * (a @ x) + >>> x0 = np.arange(12.0).reshape(6, 2) + >>> x, out = pencf(x0, obj, manifold, verbosity=0) + >>> bool(manifold.feas_eval(x) < 1e-8) + True + """ + maxit = check_maxit(maxit) + n, p = manifold._n, manifold._p + # A caller supplied starting point is used exactly as given; only + # the default one is drawn and orthonormalized. + start = ( + manifold.init_point() + if xinit is None + else as_matrix(xinit, n, p, "xinit") + ) + + x, nit, fvals, kkts, feasv, fval, kkt, fea, betout = _smopt.smpcf( + start, + _AUTO_BETA if beta is None else float(beta), + maxit, + gtol, + int(post_process), + obj_callback(obj_fun, n, p), + log_callback(verbosity, _PERIOD), + ) + out = output_dict(nit, fvals, kkts, feasv, fval, kkt, fea) + out["beta"] = betout + return x, out + + +__all__: list[str] = ["pencf"] diff --git a/src/smopt/solver/slpg.py b/src/smopt/solver/slpg.py new file mode 100644 index 0000000..87d4055 --- /dev/null +++ b/src/smopt/solver/slpg.py @@ -0,0 +1,221 @@ +"""SLPG, the penalty-free first-order solver family. + +Each entry point marshals its arguments, hands control to the Fortran 77 +driver, and turns the histories the driver fills in into the log +dictionary returned alongside the solution. +""" + +from typing import Any + +from .. import _smopt +from .._bridge import ( + Matrix, + ObjFun, + Prox, + as_matrix, + check_maxit, + log_callback, + obj_callback, + output_dict, + prox_callback, +) +from ..manifold import Stiefel + + +#: Iteration stride between the periodic progress lines. +_SMOOTH_PERIOD = 20 +_PROX_PERIOD = 50 + + +def _start(manifold: Stiefel, xinit: Matrix | None) -> Matrix: + """Resolve the starting point. + + A caller supplied one is used exactly as given; only the default is + drawn and orthonormalized. + + Args: + manifold: The manifold fixing the dimensions. + xinit: The caller's starting point, or ``None``. + + Returns: + An ``(n, p)`` matrix to start from. + """ + if xinit is None: + return manifold.init_point() + return as_matrix(xinit, manifold._n, manifold._p, "xinit") + + +def slpg_smooth( + obj_fun: ObjFun, + manifold: Stiefel, + xinit: Matrix | None = None, + maxit: int = 100, + gtol: float = 1e-5, + post_process: bool = True, # noqa: FBT001, FBT002 + verbosity: int = 2, + **kwargs: Any, # noqa: ANN401 +) -> tuple[Matrix, dict[str, object]]: + r"""Minimize a smooth objective over the Stiefel manifold. + + Args: + obj_fun: Callable mapping ``x`` to ``(fval, grad)``, where + ``grad`` is the Euclidean gradient. Returning both at once + is usually much cheaper than computing them separately. + manifold: The :class:`~smopt.manifold.Stiefel` instance fixing + the dimensions. + xinit: Starting point. A random feasible point is drawn when it + is omitted. + maxit: Maximum number of iterations. + gtol: Stationarity tolerance that stops the iteration. + post_process: Whether to round the final iterate onto the + manifold. + verbosity: ``0`` silences output, ``1`` prints the final lines, + ``2`` also prints every twentieth iteration. + **kwargs: Ignored, accepted so solvers stay interchangeable. + + Returns: + The solution and a dictionary holding the ``fvals``, ``kkts`` + and ``feas`` histories together with the final ``fval``, ``kkt`` + and ``fea`` values. + + Examples: + >>> import numpy as np + >>> from smopt import Stiefel, slpg_smooth + >>> manifold = Stiefel(6, 2) + >>> a = np.diag([5.0, 4.0, 3.0, 2.0, 1.0, 0.0]) + >>> def obj(x): + ... return float(np.sum(x * (a @ x))), 2.0 * (a @ x) + >>> x0 = np.arange(12.0).reshape(6, 2) + >>> x, out = slpg_smooth(obj, manifold, xinit=x0, verbosity=0) + >>> bool(manifold.feas_eval(x) < 1e-8) + True + """ + maxit = check_maxit(maxit) + n, p = manifold._n, manifold._p + + x, nit, fvals, kkts, feasv, fval, kkt, fea = _smopt.smslps( + _start(manifold, xinit), + maxit, + gtol, + int(post_process), + obj_callback(obj_fun, n, p), + log_callback(verbosity, _SMOOTH_PERIOD), + ) + return x, output_dict(nit, fvals, kkts, feasv, fval, kkt, fea) + + +def slpg( + obj_fun: ObjFun, + manifold: Stiefel, + xinit: Matrix | None = None, + maxit: int = 100, + prox: Prox | None = None, + gtol: float = 1e-5, + post_process: bool = True, # noqa: FBT001, FBT002 + verbosity: int = 2, + **kwargs: Any, # noqa: ANN401 +) -> tuple[Matrix, dict[str, object]]: + r"""Minimize ``f(x) + r(x)`` over the Stiefel manifold. + + The regularizer ``r`` is reached only through its proximal + operator. The multiplier of the orthogonality constraint is tracked + by an inner Arrow-Hurwicz iteration, so no penalty parameter has to + be tuned. + + Args: + obj_fun: Callable mapping ``x`` to ``(fval, grad)`` for the + smooth part ``f``. + manifold: The :class:`~smopt.manifold.Stiefel` instance fixing + the dimensions. + xinit: Starting point. A random feasible point is drawn when it + is omitted. + maxit: Maximum number of iterations. + prox: Callable mapping ``(x, eta)`` to the minimizer of + ``||y - x||_F^2 / (2 eta) + r(y)``. Defaults to the identity, + which recovers the smooth case. + gtol: Stationarity tolerance that stops the iteration. + post_process: Whether to round the final iterate onto the + manifold. + verbosity: ``0`` silences output, ``1`` prints the final lines, + ``2`` also prints every fiftieth iteration. + **kwargs: Ignored, accepted so solvers stay interchangeable. + + Returns: + The solution and a dictionary holding the ``fvals``, ``kkts`` + and ``feas`` histories together with the final ``fval``, ``kkt`` + and ``fea`` values. + """ + maxit = check_maxit(maxit) + n, p = manifold._n, manifold._p + + if prox is None: + + def prox(x: Matrix, eta: float) -> Matrix: + return x + + x, nit, fvals, kkts, feasv, fval, kkt, fea = _smopt.smslpg( + _start(manifold, xinit), + maxit, + gtol, + int(post_process), + obj_callback(obj_fun, n, p), + prox_callback(prox, n, p), + log_callback(verbosity, _PROX_PERIOD), + ) + return x, output_dict(nit, fvals, kkts, feasv, fval, kkt, fea) + + +def slpg_l21( + obj_fun: ObjFun, + manifold: Stiefel, + xinit: Matrix | None = None, + maxit: int = 100, + gamma: float = 0, + gtol: float = 1e-5, + post_process: bool = True, # noqa: FBT001, FBT002 + verbosity: int = 2, + **kwargs: Any, # noqa: ANN401 +) -> tuple[Matrix, dict[str, object]]: + r"""Minimize ``f(x) + gamma ||x||_{2,1}`` over the Stiefel manifold. + + The row-sparsity inducing :math:`\ell_{2,1}` norm has a proximal + operator and a constraint multiplier available in closed form, so + this driver needs no inner iteration. + + Args: + obj_fun: Callable mapping ``x`` to ``(fval, grad)`` for the + smooth part ``f``. + manifold: The :class:`~smopt.manifold.Stiefel` instance fixing + the dimensions. + xinit: Starting point. A random feasible point is drawn when it + is omitted. + maxit: Maximum number of iterations. + gamma: Weight of the regularization term. + gtol: Stationarity tolerance that stops the iteration. + post_process: Whether to round the final iterate onto the + manifold. + verbosity: ``0`` silences output, ``1`` prints the final lines, + ``2`` also prints every fiftieth iteration. + **kwargs: Ignored, accepted so solvers stay interchangeable. + + Returns: + The solution and a dictionary holding the ``fvals``, ``kkts`` + and ``feas`` histories together with the final ``fval``, ``kkt`` + and ``fea`` values. + """ + maxit = check_maxit(maxit) + n, p = manifold._n, manifold._p + + x, nit, fvals, kkts, feasv, fval, kkt, fea = _smopt.smsl21( + _start(manifold, xinit), + maxit, + gamma, + gtol, + int(post_process), + obj_callback(obj_fun, n, p), + log_callback(verbosity, _PROX_PERIOD), + ) + return x, output_dict(nit, fvals, kkts, feasv, fval, kkt, fea) + + +__all__: list[str] = ["slpg", "slpg_l21", "slpg_smooth"] diff --git a/src/smopt/utility/__init__.py b/src/smopt/utility/__init__.py new file mode 100644 index 0000000..66379e8 --- /dev/null +++ b/src/smopt/utility/__init__.py @@ -0,0 +1,6 @@ +"""Proximal operators of the regularizers SMOPT supports.""" + +from .utility import prox_l1, prox_l21 + + +__all__: list[str] = ["prox_l1", "prox_l21"] diff --git a/src/smopt/utility/utility.py b/src/smopt/utility/utility.py new file mode 100644 index 0000000..2088875 --- /dev/null +++ b/src/smopt/utility/utility.py @@ -0,0 +1,70 @@ +r"""Proximal operators of the regularizers SMOPT supports. + +The proximal operator of a function :math:`r` at :math:`x` with step +:math:`\eta` is the minimizer of + +.. math:: \frac{1}{2\eta} \|y - x\|_F^2 + r(y). + +Both operators below are evaluated by the Fortran 77 core, and both are +shaped so they can be handed straight to :func:`~smopt.solver.slpg`. +""" + +from .. import _smopt +from .._bridge import Matrix + + +#: Guard against dividing by a vanishing row norm. +_EPS = 1e-14 + + +def prox_l1(x: Matrix, eta: float, gamma: float = 0) -> Matrix: + r"""Proximal operator of :math:`\gamma \|x\|_1`. + + Args: + x: The point at which to evaluate the operator. + eta: The proximal step size. + gamma: Weight of the regularization term. + + Returns: + The soft-thresholded matrix, entry by entry. + + Examples: + >>> import numpy as np + >>> from smopt import prox_l1 + >>> y = prox_l1(np.array([[-3.0, 0.5, 2.0]]), 1.0, gamma=1.0) + >>> bool(np.allclose(y, [[-2.0, 0.0, 1.0]])) + True + """ + return _smopt.smpl1(x, eta, gamma) + + +def prox_l21(x: Matrix, eta: float, gamma: float = 0) -> Matrix: + r"""Proximal operator of :math:`\gamma \|x\|_{2,1}`. + + The :math:`\ell_{2,1}` norm sums the Euclidean norms of the rows of + ``x``, so the operator shrinks whole rows towards the origin and + induces row sparsity. + + Args: + x: The point at which to evaluate the operator. + eta: The proximal step size. + gamma: Weight of the regularization term. + + Returns: + The row-wise shrunk matrix. + + Examples: + >>> import numpy as np + >>> from smopt import prox_l21 + >>> x = np.array([[3.0, 4.0], [0.3, 0.4]]) + >>> bool( + ... np.allclose( + ... prox_l21(x, 1.0, gamma=1.0), [[2.4, 3.2], [0.0, 0.0]] + ... ) + ... ) + True + """ + return _smopt.smpl21(x, eta, gamma, _EPS) + + +__all__: list[str] = ["prox_l1", "prox_l21"] diff --git a/src/smpencf.f b/src/smpencf.f new file mode 100644 index 0000000..638ad76 --- /dev/null +++ b/src/smpencf.f @@ -0,0 +1,116 @@ +c----------------------------------------------------------------------- +c PENCF, a constraint dissolving penalty method for +c +c min f(X) subject to X**T*X = I_p. +c +c The search direction combines the projected gradient with a +c penalty term BETA*JC(X, C(X)) that pushes the iterate back onto +c the manifold, and feasibility is restored only once the iterate +c has drifted appreciably away from it. +c----------------------------------------------------------------------- + + SUBROUTINE SMPCF(N, P, X, BETA, MAXIT, GTOL, IPOST, OBJFUN, + $ LOGFUN, NIT, FVALS, KKTS, FEASV, FVAL, KKT, + $ FEA, BETOUT, GF, GR, GC, GRP, S, Y, XP, WP1, + $ WP2, WP3, WEIG, ROW) +c A negative BETA asks for the default 0.1*||grad f(X0)||_F, so an +c explicit BETA of zero is honoured. The value actually used is +c returned in BETOUT. + INTEGER N, P, MAXIT, IPOST, NIT + DOUBLE PRECISION X(N,P), BETA, GTOL, FVAL, KKT, FEA, BETOUT + DOUBLE PRECISION FVALS(MAXIT), KKTS(MAXIT), FEASV(MAXIT) + DOUBLE PRECISION GF(N,P), GR(N,P), GC(N,P), GRP(N,P) + DOUBLE PRECISION S(N,P), Y(N,P), XP(N,P) + DOUBLE PRECISION WP1(P,P), WP2(P,P), WP3(P,P) + DOUBLE PRECISION WEIG(P), ROW(P) + EXTERNAL OBJFUN, LOGFUN + INTEGER I, J, JJ, NP + DOUBLE PRECISION L, STEP, BT + DOUBLE PRECISION SMFRO, SMDOT, SMFEAS, SMSTEP + EXTERNAL SMFRO, SMDOT, SMFEAS, SMSTEP + + NP = N*P + NIT = 0 + KKT = 0.0D0 + FEA = 0.0D0 + + CALL OBJFUN(NP, X, FVAL, GF) + + BT = BETA + IF (BT .LT. 0.0D0) BT = 0.1D0*SMFRO(N, P, GF) + BETOUT = BT + + CALL SMJA(N, P, X, GF, GR, WP1) + CALL SMCMAP(N, P, X, WP2) + CALL SMJC(N, P, X, WP2, GC, WP1) + DO 20 J = 1, P + DO 10 I = 1, N + GR(I,J) = GR(I,J) + BT*GC(I,J) + 10 CONTINUE + 20 CONTINUE + + L = SMFRO(N, P, GF) + SMFRO(N, P, GR) + + DO 100 JJ = 1, MAXIT + IF (JJ .LE. 3) THEN + STEP = 0.01D0/L + ELSE + STEP = SMSTEP(SMDOT(N, P, S, Y), SMDOT(N, P, Y, Y), + $ 1.0D10) + END IF + + CALL SMCOPY(N, P, X, XP) + DO 40 J = 1, P + DO 30 I = 1, N + X(I,J) = X(I,J) - STEP*GR(I,J) + 30 CONTINUE + 40 CONTINUE + CALL SMFIX(N, P, X, WP1, WP2, ROW) + + DO 60 J = 1, P + DO 50 I = 1, N + S(I,J) = X(I,J) - XP(I,J) + 50 CONTINUE + 60 CONTINUE + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMCOPY(N, P, GR, GRP) + CALL SMJA(N, P, X, GF, GR, WP1) + CALL SMCMAP(N, P, X, WP2) + CALL SMJC(N, P, X, WP2, GC, WP1) + DO 80 J = 1, P + DO 70 I = 1, N + GR(I,J) = GR(I,J) + BT*GC(I,J) + Y(I,J) = GR(I,J) - GRP(I,J) + 70 CONTINUE + 80 CONTINUE + + KKT = SMFRO(N, P, GR) + FEA = SMFEAS(N, P, X, WP1) + + NIT = JJ + FVALS(JJ) = FVAL + KKTS(JJ) = KKT + FEASV(JJ) = FEA + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 0) + + IF (KKT .LT. GTOL) THEN + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 1) + GO TO 110 + END IF + 100 CONTINUE + + 110 CONTINUE + IF (IPOST .NE. 0 .AND. NIT .GE. 1) THEN + CALL SMPOST(N, P, X, WP1, WP2, WP3, WEIG, ROW) + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMJA(N, P, X, GF, GR, WP1) + KKT = SMFRO(N, P, GR) + FEA = SMFEAS(N, P, X, WP1) + CALL LOGFUN(NIT-1, FVAL, KKT, FEA, 2) + FVALS(NIT) = FVAL + KKTS(NIT) = KKT + FEASV(NIT) = FEA + END IF + RETURN + END diff --git a/src/smprox.f b/src/smprox.f new file mode 100644 index 0000000..21040ac --- /dev/null +++ b/src/smprox.f @@ -0,0 +1,84 @@ +c----------------------------------------------------------------------- +c Proximal operators of the regularizers supported by SMOPT, plus +c the closed form multiplier attached to the l_{2,1} penalty. +c +c The prox of a function r at X with step ETA is the minimizer of +c +c (1/(2*ETA))*||Y - X||_F^2 + r(Y). +c----------------------------------------------------------------------- + + SUBROUTINE SMPL1(N, P, X, ETA, GAM, Y) +c Y := prox of GAM*||.||_1 at X with step ETA, that is +c max(X - GAM*ETA, 0) + min(X + GAM*ETA, 0) entrywise. + INTEGER N, P + DOUBLE PRECISION X(N,P), Y(N,P), ETA, GAM + INTEGER I, J + DOUBLE PRECISION T + + T = GAM*ETA + DO 20 J = 1, P + DO 10 I = 1, N + Y(I,J) = DMAX1(X(I,J) - T, 0.0D0) + $ + DMIN1(X(I,J) + T, 0.0D0) + 10 CONTINUE + 20 CONTINUE + RETURN + END + + + SUBROUTINE SMPL21(N, P, X, ETA, GAM, EPS, Y) +c Y := prox of GAM*||.||_{2,1} at X with step ETA. The l_{2,1} +c norm sums the Euclidean norms of the rows of X, so the prox +c shrinks each row towards the origin. EPS guards the division by +c a vanishing row norm. + INTEGER N, P + DOUBLE PRECISION X(N,P), Y(N,P), ETA, GAM, EPS + INTEGER I, J + DOUBLE PRECISION T, R, SC + + T = GAM*ETA + DO 30 I = 1, N + R = 0.0D0 + DO 10 J = 1, P + R = R + X(I,J)*X(I,J) + 10 CONTINUE + R = DSQRT(R) + SC = DMAX1(R - T, 0.0D0)/(R + EPS) + DO 20 J = 1, P + Y(I,J) = SC*X(I,J) + 20 CONTINUE + 30 CONTINUE + RETURN + END + + + SUBROUTINE SMLM21(N, P, X, GAM, LAM) +c LAM := -GAM*X**T*diag(w)*X with w(i) = 1/(1.0D-14 + ||X(i,.)||), +c the multiplier that the l_{2,1} regularized solver attaches to +c the orthogonality constraint. + INTEGER N, P + DOUBLE PRECISION X(N,P), LAM(P,P), GAM + INTEGER I, J, K + DOUBLE PRECISION R, T + + DO 20 J = 1, P + DO 10 I = 1, P + LAM(I,J) = 0.0D0 + 10 CONTINUE + 20 CONTINUE + + DO 60 K = 1, N + R = 0.0D0 + DO 30 J = 1, P + R = R + X(K,J)*X(K,J) + 30 CONTINUE + R = 1.0D0/(1.0D-14 + DSQRT(R)) + DO 50 J = 1, P + T = GAM*R*X(K,J) + DO 40 I = 1, P + LAM(I,J) = LAM(I,J) - X(K,I)*T + 40 CONTINUE + 50 CONTINUE + 60 CONTINUE + RETURN + END diff --git a/src/smslpg.f b/src/smslpg.f new file mode 100644 index 0000000..dd7db27 --- /dev/null +++ b/src/smslpg.f @@ -0,0 +1,389 @@ +c----------------------------------------------------------------------- +c SLPG, the penalty-free first-order family of solvers for +c +c min f(X) + r(X) subject to X**T*X = I_p. +c +c Three drivers are provided: SMSLPS for a smooth f, SMSLPG for a +c general r reachable through its proximal operator, and SMSL21 for +c the l_{2,1} regularizer whose prox and multiplier are known in +c closed form. +c +c The objective is supplied by the caller through the OBJFUN +c callback, which receives the flattened iterate and returns the +c function value together with its Euclidean gradient. LOGFUN +c reports progress; the STAGE argument is 0 for an ordinary +c iteration, 1 for the line printed on convergence and 2 for the +c line printed after post-processing. +c----------------------------------------------------------------------- + + DOUBLE PRECISION FUNCTION SMSTEP(NUM, DEN, CAP) +c Barzilai-Borwein trial step |NUM/DEN| truncated at CAP. A zero +c denominator leaves the ratio unbounded, so the cap is returned. + DOUBLE PRECISION NUM, DEN, CAP + + IF (DEN .EQ. 0.0D0) THEN + SMSTEP = CAP + ELSE + SMSTEP = DMIN1(DABS(NUM/DEN), CAP) + END IF + RETURN + END + + + SUBROUTINE SMAHUR(N, P, X, G, ETA, TOL, LAM, PROXFN, + $ Z, XT, DX, WP) +c Five steps of the Arrow-Hurwicz iteration that updates the +c multiplier LAM of the orthogonality constraint at X. The loop +c stops early once the multiplier increment drops below TOL. + INTEGER N, P + DOUBLE PRECISION X(N,P), G(N,P), LAM(P,P), ETA, TOL + DOUBLE PRECISION Z(N,P), XT(N,P), DX(N,P), WP(P,P) + EXTERNAL PROXFN + INTEGER I, J, JR, NP + DOUBLE PRECISION SMFRO + EXTERNAL SMFRO + + NP = N*P + DO 20 J = 1, P + DO 10 I = 1, N + Z(I,J) = X(I,J) - ETA*G(I,J) + 10 CONTINUE + 20 CONTINUE + + DO 100 JR = 1, 5 + CALL SMJC(N, P, X, LAM, DX, WP) + DO 40 J = 1, P + DO 30 I = 1, N + XT(I,J) = Z(I,J) - ETA*DX(I,J) + 30 CONTINUE + 40 CONTINUE + CALL PROXFN(NP, XT, ETA, DX) + DO 60 J = 1, P + DO 50 I = 1, N + DX(I,J) = (DX(I,J) - X(I,J))/ETA + 50 CONTINUE + 60 CONTINUE + CALL SMJCT(N, P, X, DX, WP) + DO 80 J = 1, P + DO 70 I = 1, P + LAM(I,J) = LAM(I,J) + WP(I,J) + 70 CONTINUE + 80 CONTINUE + IF (SMFRO(P, P, WP) .LT. TOL) RETURN + 100 CONTINUE + RETURN + END + + + SUBROUTINE SMSLPS(N, P, X, MAXIT, GTOL, IPOST, OBJFUN, LOGFUN, + $ NIT, FVALS, KKTS, FEASV, FVAL, KKT, FEA, + $ GF, GR, GRP, S, Y, XP, WP1, WP2, WP3, WEIG, + $ ROW) +c SLPG for a smooth objective. A Barzilai-Borwein step is taken +c along the projected gradient and the iterate is pulled back onto +c the manifold by the feasibility restoring map A. + INTEGER N, P, MAXIT, IPOST, NIT + DOUBLE PRECISION X(N,P), GTOL, FVAL, KKT, FEA + DOUBLE PRECISION FVALS(MAXIT), KKTS(MAXIT), FEASV(MAXIT) + DOUBLE PRECISION GF(N,P), GR(N,P), GRP(N,P) + DOUBLE PRECISION S(N,P), Y(N,P), XP(N,P) + DOUBLE PRECISION WP1(P,P), WP2(P,P), WP3(P,P) + DOUBLE PRECISION WEIG(P), ROW(P) + EXTERNAL OBJFUN, LOGFUN + INTEGER I, J, JJ, NP + DOUBLE PRECISION L, STEP + DOUBLE PRECISION SMFRO, SMDOT, SMFEAS, SMSTEP + EXTERNAL SMFRO, SMDOT, SMFEAS, SMSTEP + + NP = N*P + NIT = 0 + KKT = 0.0D0 + FEA = 0.0D0 + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMJA(N, P, X, GF, GR, WP1) + L = SMFRO(N, P, GF) + SMFRO(N, P, GR) + + DO 100 JJ = 1, MAXIT + IF (JJ .LE. 3) THEN + STEP = 0.01D0/L + ELSE + STEP = SMSTEP(SMDOT(N, P, S, Y), SMDOT(N, P, Y, Y), + $ 1.0D10) + END IF + + CALL SMCOPY(N, P, X, XP) + DO 20 J = 1, P + DO 10 I = 1, N + X(I,J) = X(I,J) - STEP*GR(I,J) + 10 CONTINUE + 20 CONTINUE + CALL SMAMAP(N, P, X, WP1, WP2, ROW) + + DO 40 J = 1, P + DO 30 I = 1, N + S(I,J) = X(I,J) - XP(I,J) + 30 CONTINUE + 40 CONTINUE + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMCOPY(N, P, GR, GRP) + CALL SMJA(N, P, X, GF, GR, WP1) + DO 60 J = 1, P + DO 50 I = 1, N + Y(I,J) = GR(I,J) - GRP(I,J) + 50 CONTINUE + 60 CONTINUE + + KKT = SMFRO(N, P, GR) + FEA = SMFEAS(N, P, X, WP1) + + NIT = JJ + FVALS(JJ) = FVAL + KKTS(JJ) = KKT + FEASV(JJ) = FEA + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 0) + + IF (KKT .LT. GTOL) THEN + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 1) + GO TO 110 + END IF + 100 CONTINUE + + 110 CONTINUE + IF (IPOST .NE. 0 .AND. NIT .GE. 1) THEN + CALL SMPOST(N, P, X, WP1, WP2, WP3, WEIG, ROW) + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMJA(N, P, X, GF, GR, WP1) + KKT = SMFRO(N, P, GR) + FEA = SMFEAS(N, P, X, WP1) + CALL LOGFUN(NIT-1, FVAL, KKT, FEA, 2) + FVALS(NIT) = FVAL + KKTS(NIT) = KKT + FEASV(NIT) = FEA + END IF + RETURN + END + + + SUBROUTINE SMSLPG(N, P, X, MAXIT, GTOL, IPOST, OBJFUN, PROXFN, + $ LOGFUN, NIT, FVALS, KKTS, FEASV, FVAL, KKT, + $ FEA, GF, GR, GRAD, GRDP, S, Y, XP, Z, XT, DX, + $ LAM, WP1, WP2, WP3, WEIG, ROW, STEPS) +c SLPG for a nonsmooth regularizer reached through its proximal +c operator PROXFN. The multiplier of the orthogonality constraint +c is tracked by an Arrow-Hurwicz inner iteration so that no penalty +c parameter has to be tuned. + INTEGER N, P, MAXIT, IPOST, NIT + DOUBLE PRECISION X(N,P), GTOL, FVAL, KKT, FEA + DOUBLE PRECISION FVALS(MAXIT), KKTS(MAXIT), FEASV(MAXIT) + DOUBLE PRECISION GF(N,P), GR(N,P), GRAD(N,P), GRDP(N,P) + DOUBLE PRECISION S(N,P), Y(N,P), XP(N,P) + DOUBLE PRECISION Z(N,P), XT(N,P), DX(N,P) + DOUBLE PRECISION LAM(P,P), WP1(P,P), WP2(P,P), WP3(P,P) + DOUBLE PRECISION WEIG(P), ROW(P), STEPS(MAXIT) + EXTERNAL OBJFUN, PROXFN, LOGFUN + INTEGER I, J, K, JJ, M1, NP + DOUBLE PRECISION L, STEP, STRY, T + DOUBLE PRECISION SMFRO, SMDOT, SMFEAS, SMSTEP + EXTERNAL SMFRO, SMDOT, SMFEAS, SMSTEP + + NP = N*P + NIT = 0 + KKT = 0.0D0 + FEA = 0.0D0 + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMJA(N, P, X, GF, GR, WP1) + L = SMFRO(N, P, GF) + SMFRO(N, P, GR) + + DO 20 J = 1, P + DO 10 I = 1, P + LAM(I,J) = 0.0D0 + 10 CONTINUE + 20 CONTINUE + CALL SMAHUR(N, P, X, GR, 0.01D0/L, 0.0D0, LAM, PROXFN, + $ Z, XT, DX, WP1) + CALL SMJC(N, P, X, LAM, GRAD, WP1) + DO 40 J = 1, P + DO 30 I = 1, N + GRAD(I,J) = GR(I,J) + GRAD(I,J) + 30 CONTINUE + 40 CONTINUE + + DO 200 JJ = 1, MAXIT + IF (JJ .LE. 5) THEN + STEP = 0.01D0/L + ELSE + STEP = SMSTEP(SMDOT(N, P, S, S), SMDOT(N, P, S, Y), + $ 1.0D10) + END IF + STEPS(JJ) = STEP + + CALL SMCOPY(N, P, X, XP) + DO 60 J = 1, P + DO 50 I = 1, N + XT(I,J) = X(I,J) - STEP*GRAD(I,J) + 50 CONTINUE + 60 CONTINUE + CALL PROXFN(NP, XT, STEP, X) + CALL SMAMAP(N, P, X, WP1, WP2, ROW) + + DO 80 J = 1, P + DO 70 I = 1, N + S(I,J) = X(I,J) - XP(I,J) + 70 CONTINUE + 80 CONTINUE + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMCOPY(N, P, GRAD, GRDP) + CALL SMJA(N, P, X, GF, GR, WP1) + + M1 = MAX0(1, JJ-10) + T = 0.0D0 + DO 90 K = M1, JJ + T = T + STEPS(K) + 90 CONTINUE + STRY = T/DBLE(JJ - M1 + 1) + STRY = DMIN1(DMAX1(STRY, 1.0D-5/L), 1.0D10/L) + + FEA = SMFEAS(N, P, X, WP1) + CALL SMAHUR(N, P, X, GR, STRY, 1.0D3*FEA, LAM, PROXFN, + $ Z, XT, DX, WP1) + + CALL SMJC(N, P, X, LAM, GRAD, WP1) + DO 110 J = 1, P + DO 100 I = 1, N + GRAD(I,J) = GR(I,J) + GRAD(I,J) + Y(I,J) = GRAD(I,J) - GRDP(I,J) + 100 CONTINUE + 110 CONTINUE + + KKT = SMFRO(N, P, S)/STEP + + NIT = JJ + FVALS(JJ) = FVAL + KKTS(JJ) = KKT + FEASV(JJ) = FEA + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 0) + + IF (KKT .LT. GTOL) THEN + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 1) + GO TO 210 + END IF + 200 CONTINUE + + 210 CONTINUE + IF (IPOST .NE. 0 .AND. NIT .GE. 1) THEN + CALL SMPOST(N, P, X, WP1, WP2, WP3, WEIG, ROW) + CALL OBJFUN(NP, X, FVAL, GF) + FEA = SMFEAS(N, P, X, WP1) + CALL LOGFUN(NIT-1, FVAL, KKT, FEA, 2) + FVALS(NIT) = FVAL + KKTS(NIT) = KKT + FEASV(NIT) = FEA + END IF + RETURN + END + + + SUBROUTINE SMSL21(N, P, X, MAXIT, GAM, GTOL, IPOST, OBJFUN, + $ LOGFUN, NIT, FVALS, KKTS, FEASV, FVAL, KKT, + $ FEA, GF, GR, GRAD, GRDP, S, Y, XP, XT, LAM, + $ WP1, WP2, WP3, WEIG, ROW) +c SLPG for the l_{2,1} regularized objective f(X) + GAM*||X||_{2,1}. +c Both the prox and the constraint multiplier are available in +c closed form, so no inner iteration is needed. + INTEGER N, P, MAXIT, IPOST, NIT + DOUBLE PRECISION X(N,P), GAM, GTOL, FVAL, KKT, FEA + DOUBLE PRECISION FVALS(MAXIT), KKTS(MAXIT), FEASV(MAXIT) + DOUBLE PRECISION GF(N,P), GR(N,P), GRAD(N,P), GRDP(N,P) + DOUBLE PRECISION S(N,P), Y(N,P), XP(N,P), XT(N,P) + DOUBLE PRECISION LAM(P,P), WP1(P,P), WP2(P,P), WP3(P,P) + DOUBLE PRECISION WEIG(P), ROW(P) + EXTERNAL OBJFUN, LOGFUN + INTEGER I, J, JJ, NP + DOUBLE PRECISION L, STEP + DOUBLE PRECISION SMFRO, SMDOT, SMFEAS, SMSTEP + EXTERNAL SMFRO, SMDOT, SMFEAS, SMSTEP + + NP = N*P + NIT = 0 + KKT = 0.0D0 + FEA = 0.0D0 + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMJA(N, P, X, GF, GR, WP1) + L = SMFRO(N, P, GF) + SMFRO(N, P, GR) + + CALL SMLM21(N, P, X, GAM, LAM) + CALL SMJC(N, P, X, LAM, GRAD, WP1) + DO 20 J = 1, P + DO 10 I = 1, N + GRAD(I,J) = GR(I,J) + GRAD(I,J) + 10 CONTINUE + 20 CONTINUE + + DO 100 JJ = 1, MAXIT + IF (JJ .LE. 5) THEN + STEP = 0.001D0/L + ELSE + STEP = SMSTEP(SMDOT(N, P, S, S), SMDOT(N, P, S, Y), + $ 1.0D5) + END IF + + CALL SMCOPY(N, P, X, XP) + DO 40 J = 1, P + DO 30 I = 1, N + XT(I,J) = X(I,J) - STEP*GRAD(I,J) + 30 CONTINUE + 40 CONTINUE + CALL SMPL21(N, P, XT, STEP, GAM, 1.0D-16, X) + CALL SMAMAP(N, P, X, WP1, WP2, ROW) + + DO 60 J = 1, P + DO 50 I = 1, N + S(I,J) = X(I,J) - XP(I,J) + 50 CONTINUE + 60 CONTINUE + + CALL OBJFUN(NP, X, FVAL, GF) + CALL SMCOPY(N, P, GRAD, GRDP) + CALL SMJA(N, P, X, GF, GR, WP1) + + CALL SMLM21(N, P, X, GAM, LAM) + CALL SMJC(N, P, X, LAM, GRAD, WP1) + DO 80 J = 1, P + DO 70 I = 1, N + GRAD(I,J) = GR(I,J) + GRAD(I,J) + Y(I,J) = GRAD(I,J) - GRDP(I,J) + 70 CONTINUE + 80 CONTINUE + + KKT = SMFRO(N, P, S)/STEP + FEA = SMFEAS(N, P, X, WP1) + + NIT = JJ + FVALS(JJ) = FVAL + KKTS(JJ) = KKT + FEASV(JJ) = FEA + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 0) + + IF (KKT .LT. GTOL) THEN + CALL LOGFUN(JJ-1, FVAL, KKT, FEA, 1) + GO TO 110 + END IF + 100 CONTINUE + + 110 CONTINUE + IF (IPOST .NE. 0 .AND. NIT .GE. 1) THEN + CALL SMPOST(N, P, X, WP1, WP2, WP3, WEIG, ROW) + CALL OBJFUN(NP, X, FVAL, GF) + FEA = SMFEAS(N, P, X, WP1) + CALL LOGFUN(NIT-1, FVAL, KKT, FEA, 2) + FVALS(NIT) = FVAL + KKTS(NIT) = KKT + FEASV(NIT) = FEA + END IF + RETURN + END diff --git a/tests/.gitkeep b/tests/.gitkeep deleted file mode 100644 index e69de29..0000000 diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 0000000..149609e --- /dev/null +++ b/tests/conftest.py @@ -0,0 +1,35 @@ +"""Shared fixtures and problem builders for the SMOPT test suite.""" + +import numpy as np +import pytest +import reference + +import smopt + + +@pytest.fixture +def rng() -> np.random.Generator: + """A seeded generator so every test is reproducible.""" + return np.random.default_rng(20260825) + + +def orthonormal(rng: np.random.Generator, n: int, p: int) -> np.ndarray: + """Draw a random point on the Stiefel manifold.""" + q, _ = np.linalg.qr(rng.standard_normal((n, p))) + return np.ascontiguousarray(q[:, :p]) + + +def quadratic(a_diag: np.ndarray): + """Build ``tr(X^T A X)`` and its gradient for a diagonal ``A``.""" + a = np.asarray(a_diag, dtype=np.float64) + + def obj(x: np.ndarray) -> tuple[float, np.ndarray]: + ax = a[:, None] * x + return float(np.sum(x * ax)), 2.0 * ax + + return obj + + +def manifolds(n: int, p: int) -> tuple[smopt.Stiefel, reference.Stiefel]: + """Return the shipped manifold and its reference twin.""" + return smopt.Stiefel(n, p), reference.Stiefel(n, p) diff --git a/tests/reference.py b/tests/reference.py new file mode 100644 index 0000000..76ef608 --- /dev/null +++ b/tests/reference.py @@ -0,0 +1,411 @@ +"""A NumPy transcription of the algorithm SMOPT ports to Fortran 77. + +This module exists purely as a test oracle. It mirrors the original +implementation statement by statement so that the Fortran drivers can be +checked against it, and it is deliberately kept free of the validation, +packaging and reporting concerns that the shipped package handles. + +Names follow the conventions of this project rather than those of the +original, so the oracle and the package under test can be driven through +the same attribute lookups. The one behavioural departure is +:meth:`Stiefel.init_point`, whose ``xinit == None`` test raises on array +input; it is spelled ``is None`` here so the reference can be driven from +the tests. +""" + +import numpy as np +from numpy.linalg import norm, svd + + +class Stiefel: + """Reference geometry of the Stiefel manifold.""" + + def __init__(self, n, p): + self._n = n + self._p = p + self.dim = n * p + + def phi(self, m): + return (m + m.T) / 2 + + def a(self, x): + xx = x.T @ x + feas_tmp = norm(xx - np.eye(self._p), "fro") + if feas_tmp < 0.5: + return 1.5 * x - x @ (xx / 2) + return np.linalg.solve((xx + np.eye(self._p)) / 2, x.T).T + + def ja(self, x, g): + return g - x @ self.phi(x.T @ g) + + def jc(self, x, lam): + return x @ self.phi(lam) + + def jc_transpose(self, x, d): + return self.phi(x.T @ d) + + def c(self, x): + return x.T @ x - np.eye(self._p) + + def feas_eval(self, x): + return norm(self.c(x), "fro") + + def init_point(self, xinit=None): + if xinit is None: + xinit = np.random.randn(self._n, self._p) + if norm(xinit.T @ xinit - np.eye(self._p), "fro") > 1e-6: + xinit, _ = np.linalg.qr(xinit) + return xinit + + def post_process(self, x): + ux, _, vx = svd(x, full_matrices=False) + return ux @ vx + + +def slpg_smooth( + obj_fun, + manifold, + xinit=None, + maxit=100, + gtol=1e-5, + post_process=True, + verbosity=0, +): + """Reference SLPG for a smooth objective.""" + kkts, feas, fvals = [], [], [] + + if xinit is None: + xinit = manifold.init_point() + + x = xinit + fval, gradf = obj_fun(x) + gradr = manifold.ja(x, gradf) + lip = norm(gradf, "fro") + norm(gradr, "fro") + + s = y = None + for jj in range(maxit): + if jj < 3: + stepsize = 0.01 / lip + else: + stepsize = np.abs(np.sum(s * y) / np.sum(y * y)) + stepsize = np.min((stepsize, 1e10)) + + x_p = x + x = x - stepsize * gradr + x = manifold.a(x) + s = x - x_p + + fval, gradf = obj_fun(x) + gradr_p = gradr + gradr = manifold.ja(x, gradf) + y = gradr - gradr_p + + substationarity = norm(gradr, "fro") + feasibility = manifold.feas_eval(x) + + kkts.append(substationarity) + feas.append(feasibility) + fvals.append(fval) + + if substationarity < gtol: + break + + if post_process: + x = manifold.post_process(x) + fval, gradf = obj_fun(x) + gradr = manifold.ja(x, gradf) + substationarity = norm(gradr, "fro") + feasibility = manifold.feas_eval(x) + kkts[-1] = substationarity + feas[-1] = feasibility + fvals[-1] = fval + + return x, { + "kkts": kkts, + "fvals": fvals, + "fea": feasibility, + "kkt": substationarity, + "fval": fval, + "feas": feas, + } + + +def arrow_hurwicz_slpg(x, g, eta, prox, lam, manifold, tol=0): + """Reference Arrow-Hurwicz multiplier update.""" + lam_temp = lam + try_stepsize = eta + z_tmp = x - try_stepsize * g + for _ in range(5): + x_try = prox( + z_tmp - try_stepsize * manifold.jc(x, lam_temp), try_stepsize + ) + d_x = 1 / try_stepsize * (x_try - x) + lam_inc = manifold.jc_transpose(x, d_x) + lam_temp = lam_temp + lam_inc + if norm(lam_inc, "fro") < tol: + break + return lam_temp + + +def slpg( + obj_fun, + manifold, + xinit=None, + maxit=100, + prox=lambda x, eta: x, + gtol=1e-5, + post_process=True, + verbosity=0, +): + """Reference SLPG for a proximable regularizer.""" + kkts, feas, fvals, steps = [], [], [], [] + + if xinit is None: + xinit = manifold.init_point() + + p = manifold._p + x = xinit + fval, gradf = obj_fun(x) + gradr = manifold.ja(x, gradf) + lip = norm(gradf, "fro") + norm(gradr, "fro") + + lam = np.zeros([p, p]) + lam = arrow_hurwicz_slpg(x, gradr, 0.01 / lip, prox, lam, manifold) + grad = gradr + manifold.jc(x, lam) + + s = y = None + for jj in range(maxit): + if jj < 5: + stepsize = 0.01 / lip + else: + stepsize = np.abs(np.sum(s * s) / np.sum(s * y)) + stepsize = np.min((stepsize, 1e10)) + + x_p = x + steps.append(stepsize) + + x = prox(x - stepsize * (gradr + manifold.jc(x, lam)), stepsize) + x = manifold.a(x) + s = x - x_p + + fval, gradf = obj_fun(x) + grad_p = grad + gradr = manifold.ja(x, gradf) + + stepsize_try = np.average(steps[np.maximum(0, jj - 10) :]) + stepsize_try = np.minimum( + np.maximum(stepsize_try, 1e-5 / lip), 1e10 / lip + ) + + tol_aw = 1000 * manifold.feas_eval(x) + lam = arrow_hurwicz_slpg( + x, gradr, stepsize_try, prox, lam, manifold, tol=tol_aw + ) + grad = gradr + manifold.jc(x, lam) + y = grad - grad_p + + substationarity = norm(s / stepsize, "fro") + feasibility = manifold.feas_eval(x) + + kkts.append(substationarity) + feas.append(feasibility) + fvals.append(fval) + + if substationarity < gtol: + break + + if post_process: + x = manifold.post_process(x) + fval, gradf = obj_fun(x) + feasibility = manifold.feas_eval(x) + kkts[-1] = substationarity + feas[-1] = feasibility + fvals[-1] = fval + + return x, { + "kkts": kkts, + "fvals": fvals, + "fea": feasibility, + "kkt": substationarity, + "fval": fval, + "feas": feas, + } + + +def slpg_l21( + obj_fun, + manifold, + xinit=None, + maxit=100, + gamma=0, + gtol=1e-5, + post_process=True, + verbosity=0, +): + """Reference SLPG for the l_{2,1} regularized objective.""" + + def prox(x_input, eta): + x_ref = np.sqrt(np.sum(x_input**2, axis=1, keepdims=True)) + x_ref_reduce = np.maximum(x_ref - gamma * eta, 0) + return (x_ref_reduce / (x_ref + 1e-16)) * x_input + + def generate_lam(x_input): + x_ref = 1 / (1e-14 + np.sqrt(np.sum(x_input**2, axis=1, keepdims=True))) + return -x_input.T @ (x_ref * x_input) + + kkts, feas, fvals = [], [], [] + + if xinit is None: + xinit = manifold.init_point() + + x = xinit + fval, gradf = obj_fun(x) + gradr = manifold.ja(x, gradf) + lip = norm(gradf, "fro") + norm(gradr, "fro") + + lam = gamma * generate_lam(x) + grad = gradr + manifold.jc(x, lam) + + s = y = None + for jj in range(maxit): + if jj < 5: + stepsize = 0.001 / lip + else: + stepsize = np.abs(np.sum(s * s) / np.sum(s * y)) + stepsize = np.min((stepsize, 1e5)) + + x_p = x + x = prox(x - stepsize * grad, stepsize) + x = manifold.a(x) + s = x - x_p + + fval, gradf = obj_fun(x) + grad_p = grad + gradr = manifold.ja(x, gradf) + + lam = gamma * generate_lam(x) + grad = gradr + manifold.jc(x, lam) + y = grad - grad_p + + substationarity = norm(s / stepsize, "fro") + feasibility = manifold.feas_eval(x) + + kkts.append(substationarity) + feas.append(feasibility) + fvals.append(fval) + + if substationarity < gtol: + break + + if post_process: + x = manifold.post_process(x) + fval, gradf = obj_fun(x) + feasibility = manifold.feas_eval(x) + kkts[-1] = substationarity + feas[-1] = feasibility + fvals[-1] = fval + + return x, { + "kkts": kkts, + "fvals": fvals, + "fea": feasibility, + "kkt": substationarity, + "fval": fval, + "feas": feas, + } + + +def pencf( + xinit, + obj_fun, + manifold, + beta=None, + maxit=100, + gtol=1e-5, + post_process=True, + verbosity=0, +): + """Reference PenCF.""" + kkts, feas, fvals = [], [], [] + + p = manifold._p + x = xinit + fval, gradf = obj_fun(x) + + if beta is None: + beta = 0.1 * norm(gradf, "fro") + + gradr = manifold.ja(x, gradf) + beta * manifold.jc(x, manifold.c(x)) + lip = norm(gradf, "fro") + norm(gradr, "fro") + + s = y = None + for jj in range(maxit): + if jj < 3: + stepsize = 0.01 / lip + else: + stepsize = np.abs(np.sum(s * y) / np.sum(y * y)) + stepsize = np.min((stepsize, 1e10)) + + x_p = x + x = x - stepsize * gradr + + xx = x.T @ x + feas_tmp = manifold.feas_eval(x) + if feas_tmp > 1e-1: + if feas_tmp < 0.5: + x = 1.5 * x - x @ (xx / 2) + else: + x = np.linalg.solve((xx + np.eye(p)) / 2, x.T).T + + if norm(x, "fro") > 1.001 * np.sqrt(p): + x = x * (1.001 * np.sqrt(p) / norm(x, "fro")) + + s = x - x_p + + fval, gradf = obj_fun(x) + gradr_p = gradr + gradr = manifold.ja(x, gradf) + beta * manifold.jc(x, manifold.c(x)) + y = gradr - gradr_p + + substationarity = norm(gradr, "fro") + feasibility = manifold.feas_eval(x) + + kkts.append(substationarity) + feas.append(feasibility) + fvals.append(fval) + + if substationarity < gtol: + break + + if post_process: + x = manifold.post_process(x) + fval, gradf = obj_fun(x) + gradr = manifold.ja(x, gradf) + substationarity = norm(gradr, "fro") + feasibility = manifold.feas_eval(x) + kkts[-1] = substationarity + feas[-1] = feasibility + fvals[-1] = fval + + return x, { + "kkts": kkts, + "fvals": fvals, + "fea": feasibility, + "kkt": substationarity, + "fval": fval, + "feas": feas, + } + + +def prox_l1(x_input, eta, gamma=0): + """Reference proximal operator of the l_1 norm.""" + return np.maximum(x_input - gamma * eta, 0) + np.minimum( + x_input + gamma * eta, 0 + ) + + +def prox_l21(x_input, eta, gamma=0): + """Reference proximal operator of the l_{2,1} norm.""" + x_ref = np.sqrt(np.sum(x_input**2, axis=1, keepdims=True)) + x_ref_reduce = np.maximum(x_ref - gamma * eta, 0) + return x_ref_reduce / (x_ref + 1e-14) * x_input diff --git a/tests/test_manifold.py b/tests/test_manifold.py new file mode 100644 index 0000000..49f9143 --- /dev/null +++ b/tests/test_manifold.py @@ -0,0 +1,172 @@ +"""The Fortran manifold maps must agree with the NumPy reference.""" + +import numpy as np +import pytest +import reference +from conftest import orthonormal + +import smopt + + +SHAPES = [(1, 1), (5, 1), (12, 3), (40, 7), (9, 9)] + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_phi_symmetrizes(rng: np.random.Generator, n: int, p: int) -> None: + m = smopt.Stiefel(n, p) + a = rng.standard_normal((p, p)) + got = m.phi(a) + assert np.allclose(got, reference.Stiefel(n, p).phi(a)) + assert np.allclose(got, got.T) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_phi_does_not_mutate_input( + rng: np.random.Generator, n: int, p: int +) -> None: + m = smopt.Stiefel(n, p) + a = rng.standard_normal((p, p)) + before = a.copy() + m.phi(a) + assert np.array_equal(a, before) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_constraint_and_feasibility( + rng: np.random.Generator, n: int, p: int +) -> None: + m, ref = smopt.Stiefel(n, p), reference.Stiefel(n, p) + x = rng.standard_normal((n, p)) + assert np.allclose(m.c(x), ref.c(x)) + assert m.feas_eval(x) == pytest.approx(ref.feas_eval(x)) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_feasibility_vanishes_on_the_manifold( + rng: np.random.Generator, n: int, p: int +) -> None: + m = smopt.Stiefel(n, p) + assert m.feas_eval(orthonormal(rng, n, p)) < 1e-12 + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_jacobians(rng: np.random.Generator, n: int, p: int) -> None: + m, ref = smopt.Stiefel(n, p), reference.Stiefel(n, p) + x = orthonormal(rng, n, p) + g = rng.standard_normal((n, p)) + d = rng.standard_normal((n, p)) + lam = rng.standard_normal((p, p)) + + assert np.allclose(m.ja(x, g), ref.ja(x, g)) + assert np.allclose(m.jc(x, lam), ref.jc(x, lam)) + assert np.allclose(m.jc_transpose(x, d), ref.jc_transpose(x, d)) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_ja_output_is_tangent(rng: np.random.Generator, n: int, p: int) -> None: + """``JA`` removes the symmetric part of ``X^T G``.""" + m = smopt.Stiefel(n, p) + x = orthonormal(rng, n, p) + g = rng.standard_normal((n, p)) + r = m.ja(x, g) + xtr = x.T @ r + assert np.allclose(xtr, -xtr.T, atol=1e-10) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +@pytest.mark.parametrize("scale", [1e-3, 0.05, 0.4, 3.0]) +def test_a_map_matches_reference( + rng: np.random.Generator, n: int, p: int, scale: float +) -> None: + """Both the near-manifold expansion and the exact solve branch.""" + m, ref = smopt.Stiefel(n, p), reference.Stiefel(n, p) + x = orthonormal(rng, n, p) + scale * rng.standard_normal((n, p)) + assert np.allclose(m.a(x), ref.a(x), atol=1e-10) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_a_map_improves_feasibility( + rng: np.random.Generator, n: int, p: int +) -> None: + m = smopt.Stiefel(n, p) + x = orthonormal(rng, n, p) + 0.05 * rng.standard_normal((n, p)) + assert m.feas_eval(m.a(x)) < m.feas_eval(x) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_post_process_is_the_polar_factor( + rng: np.random.Generator, n: int, p: int +) -> None: + m, ref = smopt.Stiefel(n, p), reference.Stiefel(n, p) + x = rng.standard_normal((n, p)) + got = m.post_process(x) + assert np.allclose(got, ref.post_process(x), atol=1e-9) + assert m.feas_eval(got) < 1e-12 + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_init_point_lands_on_the_manifold( + rng: np.random.Generator, n: int, p: int +) -> None: + m = smopt.Stiefel(n, p) + assert m.feas_eval(m.init_point(rng.standard_normal((n, p)))) < 1e-12 + assert m.feas_eval(m.init_point()) < 1e-12 + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +def test_init_point_keeps_a_feasible_argument( + rng: np.random.Generator, n: int, p: int +) -> None: + m = smopt.Stiefel(n, p) + x = orthonormal(rng, n, p) + assert np.allclose(m.init_point(x), x) + + +def test_rejects_bad_dimensions() -> None: + with pytest.raises(ValueError, match="must be positive"): + smopt.Stiefel(0, 1) + with pytest.raises(ValueError, match="must not exceed"): + smopt.Stiefel(2, 3) + + +def test_rejects_mismatched_shapes(rng: np.random.Generator) -> None: + m = smopt.Stiefel(6, 2) + with pytest.raises(ValueError, match="must have shape"): + m.c(rng.standard_normal((6, 3))) + with pytest.raises(ValueError, match="must have shape"): + m.jc(orthonormal(rng, 6, 2), rng.standard_normal((3, 3))) + + +def test_dim_attribute() -> None: + assert smopt.Stiefel(7, 3).dim == 21 + + +@pytest.mark.parametrize("rank", [0, 1, 2]) +def test_post_process_completes_a_rank_deficient_point( + rng: np.random.Generator, rank: int +) -> None: + """A vanishing singular value leaves its direction unconstrained. + + There is no polar factor then, so the missing directions get an + arbitrary orthonormal completion rather than dividing by zero. + """ + n, p = 9, 3 + m = smopt.Stiefel(n, p) + x = np.zeros((n, p)) + if rank: + x[:, :rank] = orthonormal(rng, n, rank) + + got = m.post_process(x) + + assert np.all(np.isfinite(got)) + assert m.feas_eval(got) < 1e-12 + # Whatever the completion picks, it must not disturb the directions + # the input did pin down. + assert np.allclose(got @ got.T @ x, x, atol=1e-10) + + +def test_post_process_of_a_zero_matrix_is_still_feasible() -> None: + m = smopt.Stiefel(5, 2) + got = m.post_process(np.zeros((5, 2))) + assert np.all(np.isfinite(got)) + assert m.feas_eval(got) < 1e-12 diff --git a/tests/test_solvers.py b/tests/test_solvers.py new file mode 100644 index 0000000..bfbcd90 --- /dev/null +++ b/tests/test_solvers.py @@ -0,0 +1,344 @@ +"""The Fortran solver drivers must reproduce the reference iteration. + +Each solver is run side by side with the NumPy transcription in +``reference.py`` from the same starting point, and the whole trajectory +is compared, not just the final answer. +""" + +import numpy as np +import pytest +import reference +from conftest import orthonormal, quadratic + +import smopt + + +# The two implementations are transcriptions of one another, not +# bit-identical: NumPy sums through BLAS while the Fortran kernels sum in +# their own order. That seed then amplifies, because the +# Barzilai-Borwein step divides differences of nearly equal quantities +# and the l_{2,1} prox thresholds whole rows on or off. Iterates agree to +# ~1e-16 at the start, ~6e-10 by iteration 10 and only ~3e-7 by iteration +# 15. These tests therefore compare a short horizon strictly; that the +# solvers actually converge is checked separately below. +TRACE_MAXIT = 10 +TRACE_RTOL = 1e-7 +TRACE_ATOL = 1e-10 +TRACE_X_ATOL = 1e-7 + + +def eig_problem(rng: np.random.Generator, n: int, p: int): + """A trace minimization problem with a known optimal value.""" + a = np.sort(rng.uniform(0.5, 10.0, size=n)) + return quadratic(a), float(np.sum(a[:p])) + + +def assert_same_trace(got: dict, want: dict) -> None: + """Compare a Fortran trajectory with the reference one.""" + kw = {"rtol": TRACE_RTOL, "atol": TRACE_ATOL} + assert len(got["fvals"]) == len(want["fvals"]) + assert np.allclose(got["fvals"], want["fvals"], **kw) + assert np.allclose(got["kkts"], want["kkts"], **kw) + assert np.allclose(got["feas"], want["feas"], **kw) + scalar = {"rel": TRACE_RTOL, "abs": TRACE_ATOL} + assert got["fval"] == pytest.approx(want["fval"], **scalar) + assert got["kkt"] == pytest.approx(want["kkt"], **scalar) + assert got["fea"] == pytest.approx(want["fea"], **scalar) + + +@pytest.mark.parametrize(("n", "p"), [(20, 3), (60, 5), (8, 8)]) +@pytest.mark.parametrize("post_process", [True, False]) +def test_slpg_smooth_matches_reference( + rng: np.random.Generator, n: int, p: int, post_process: bool +) -> None: + m, ref_m = smopt.Stiefel(n, p), reference.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + got_x, got = smopt.slpg_smooth( + obj, + m, + xinit=x0.copy(), + maxit=TRACE_MAXIT, + post_process=post_process, + verbosity=0, + ) + want_x, want = reference.slpg_smooth( + obj, + ref_m, + xinit=x0.copy(), + maxit=TRACE_MAXIT, + post_process=post_process, + ) + + assert_same_trace(got, want) + assert np.allclose(got_x, want_x, atol=TRACE_X_ATOL) + + +@pytest.mark.parametrize(("n", "p"), [(20, 3), (60, 5)]) +def test_slpg_matches_reference( + rng: np.random.Generator, n: int, p: int +) -> None: + m, ref_m = smopt.Stiefel(n, p), reference.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + gamma = 0.05 + + def prox(x: np.ndarray, eta: float) -> np.ndarray: + return smopt.prox_l1(x, eta, gamma=gamma) + + def ref_prox(x: np.ndarray, eta: float) -> np.ndarray: + return reference.prox_l1(x, eta, gamma=gamma) + + got_x, got = smopt.slpg( + obj, m, xinit=x0.copy(), maxit=TRACE_MAXIT, prox=prox, verbosity=0 + ) + want_x, want = reference.slpg( + obj, ref_m, xinit=x0.copy(), maxit=TRACE_MAXIT, prox=ref_prox + ) + + assert_same_trace(got, want) + assert np.allclose(got_x, want_x, atol=TRACE_X_ATOL) + + +@pytest.mark.parametrize(("n", "p"), [(20, 3), (60, 5)]) +def test_slpg_without_a_prox_matches_reference( + rng: np.random.Generator, n: int, p: int +) -> None: + """The default prox is the identity, recovering the smooth case.""" + m, ref_m = smopt.Stiefel(n, p), reference.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + got_x, got = smopt.slpg( + obj, m, xinit=x0.copy(), maxit=TRACE_MAXIT, verbosity=0 + ) + want_x, want = reference.slpg( + obj, ref_m, xinit=x0.copy(), maxit=TRACE_MAXIT + ) + + assert_same_trace(got, want) + assert np.allclose(got_x, want_x, atol=TRACE_X_ATOL) + + +@pytest.mark.parametrize(("n", "p"), [(20, 3), (60, 5)]) +@pytest.mark.parametrize("gamma", [0.0, 0.02]) +def test_slpg_l21_matches_reference( + rng: np.random.Generator, n: int, p: int, gamma: float +) -> None: + m, ref_m = smopt.Stiefel(n, p), reference.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + got_x, got = smopt.slpg_l21( + obj, m, xinit=x0.copy(), maxit=TRACE_MAXIT, gamma=gamma, verbosity=0 + ) + want_x, want = reference.slpg_l21( + obj, ref_m, xinit=x0.copy(), maxit=TRACE_MAXIT, gamma=gamma + ) + + assert_same_trace(got, want) + assert np.allclose(got_x, want_x, atol=TRACE_X_ATOL) + + +@pytest.mark.parametrize(("n", "p"), [(20, 3), (60, 5)]) +@pytest.mark.parametrize("beta", [None, 0.5]) +def test_pencf_matches_reference( + rng: np.random.Generator, n: int, p: int, beta: float | None +) -> None: + m, ref_m = smopt.Stiefel(n, p), reference.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + got_x, got = smopt.pencf( + x0.copy(), obj, m, beta=beta, maxit=TRACE_MAXIT, verbosity=0 + ) + want_x, want = reference.pencf( + x0.copy(), obj, ref_m, beta=beta, maxit=TRACE_MAXIT + ) + + assert_same_trace(got, want) + assert np.allclose(got_x, want_x, atol=TRACE_X_ATOL) + + +def test_pencf_reports_the_beta_it_used(rng: np.random.Generator) -> None: + n, p = 20, 3 + m = smopt.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + _, out = smopt.pencf(x0.copy(), obj, m, maxit=5, verbosity=0) + _, grad = obj(x0) + assert out["beta"] == pytest.approx(0.1 * np.linalg.norm(grad, "fro")) + + _, out = smopt.pencf(x0.copy(), obj, m, beta=0.0, maxit=5, verbosity=0) + assert out["beta"] == pytest.approx(0.0) + + +SOLVERS = ["slpg_smooth", "slpg", "slpg_l21"] + + +@pytest.mark.parametrize("name", SOLVERS) +def test_solvers_reach_the_known_minimum( + rng: np.random.Generator, name: str +) -> None: + """Trace minimization converges to the sum of the smallest values.""" + n, p = 50, 4 + m = smopt.Stiefel(n, p) + obj, best = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + solver = getattr(smopt, name) + x, out = solver(obj, m, xinit=x0, maxit=3000, gtol=1e-9, verbosity=0) + + assert m.feas_eval(x) < 1e-10 + assert out["fval"] == pytest.approx(best, rel=1e-5) + + +#: pencf lands near the optimum rather than on it. Holding the problem +#: below fixed and varying only the starting point over 60 draws, the +#: relative error has median 7.6e-07 and worst case 1.2e-04; a 150-seed +#: sweep over problems too tops out at 1.0e-04, and the NumPy reference +#: in tests/reference.py is worse still at 1.7e-04. That is the penalty +#: method, not this port. The slpg solvers above reach machine precision +#: on the same instance, hence their much tighter tolerance. +PENALTY_RTOL = 1e-3 + + +def test_pencf_reaches_the_known_minimum(rng: np.random.Generator) -> None: + """pencf converges to the optimum to within a penalty method's reach.""" + n, p = 50, 4 + m = smopt.Stiefel(n, p) + obj, best = eig_problem(rng, n, p) + + x, out = smopt.pencf( + orthonormal(rng, n, p), obj, m, maxit=3000, gtol=1e-9, verbosity=0 + ) + + # Feasibility is the part the method does guarantee: it holds to + # ~3e-15 across every instance measured, so it is asserted tightly. + assert m.feas_eval(x) < 1e-10 + assert out["fval"] == pytest.approx(best, rel=PENALTY_RTOL) + # A feasible point can never beat the true minimum; catching that + # would mean the objective and the constraint had come apart. + assert out["fval"] > best - 1e-8 + + +def test_l21_regularization_induces_row_sparsity( + rng: np.random.Generator, +) -> None: + n, p = 40, 3 + m = smopt.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + dense, _ = smopt.slpg_l21( + obj, m, xinit=x0.copy(), maxit=500, gamma=0.0, verbosity=0 + ) + sparse, _ = smopt.slpg_l21( + obj, m, xinit=x0.copy(), maxit=500, gamma=1.0, verbosity=0 + ) + + def live_rows(x: np.ndarray) -> int: + return int(np.sum(np.linalg.norm(x, axis=1) > 1e-6)) + + assert live_rows(sparse) < live_rows(dense) + + +@pytest.mark.parametrize("name", SOLVERS) +def test_histories_have_one_entry_per_iteration( + rng: np.random.Generator, name: str +) -> None: + n, p = 15, 2 + m = smopt.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + + solver = getattr(smopt, name) + _, out = solver( + obj, m, xinit=orthonormal(rng, n, p), maxit=7, gtol=0.0, verbosity=0 + ) + + assert len(out["fvals"]) == 7 + assert len(out["kkts"]) == 7 + assert len(out["feas"]) == 7 + + +@pytest.mark.parametrize("name", SOLVERS) +def test_a_random_start_is_drawn_when_none_is_given( + rng: np.random.Generator, name: str +) -> None: + n, p = 15, 2 + m = smopt.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + + x, _ = getattr(smopt, name)(obj, m, maxit=20, verbosity=0) + assert x.shape == (n, p) + assert m.feas_eval(x) < 1e-10 + + +@pytest.mark.parametrize("name", SOLVERS) +def test_rejects_an_empty_iteration_budget( + rng: np.random.Generator, name: str +) -> None: + m = smopt.Stiefel(6, 2) + obj, _ = eig_problem(rng, 6, 2) + with pytest.raises(ValueError, match="maxit must be at least 1"): + getattr(smopt, name)(obj, m, maxit=0, verbosity=0) + + +def test_verbosity_controls_printing( + rng: np.random.Generator, capsys: pytest.CaptureFixture +) -> None: + n, p = 15, 2 + m = smopt.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + x0 = orthonormal(rng, n, p) + + smopt.slpg_smooth(obj, m, xinit=x0.copy(), maxit=25, verbosity=0) + assert capsys.readouterr().out == "" + + smopt.slpg_smooth(obj, m, xinit=x0.copy(), maxit=25, verbosity=2) + out = capsys.readouterr().out + assert "Iter:0" in out + assert "Iter:20" in out + assert "Post-processing" in out + + +def test_the_objective_sees_the_iterate_as_a_matrix( + rng: np.random.Generator, +) -> None: + """The callback must receive an ``(n, p)`` array, not a flat one.""" + n, p = 12, 3 + m = smopt.Stiefel(n, p) + seen = [] + + def obj(x: np.ndarray) -> tuple[float, np.ndarray]: + seen.append(x.shape) + return float(np.sum(x * x)), 2.0 * x + + smopt.slpg_smooth( + obj, m, xinit=orthonormal(rng, n, p), maxit=3, verbosity=0 + ) + assert seen + assert set(seen) == {(n, p)} + + +def test_l21_stays_feasible_when_it_collapses_the_rank( + rng: np.random.Generator, +) -> None: + """A heavy penalty can zero enough rows to make the iterate singular. + + Post-processing must still return a point on the manifold instead of + dividing by a vanishing singular value. + """ + n, p = 20, 3 + m = smopt.Stiefel(n, p) + obj, _ = eig_problem(rng, n, p) + + x, out = smopt.slpg_l21( + obj, m, xinit=orthonormal(rng, n, p), maxit=40, gamma=0.1, verbosity=0 + ) + + assert np.all(np.isfinite(x)) + assert m.feas_eval(x) < 1e-10 + assert np.isfinite(out["fval"]) diff --git a/tests/test_utility.py b/tests/test_utility.py new file mode 100644 index 0000000..b54294a --- /dev/null +++ b/tests/test_utility.py @@ -0,0 +1,77 @@ +"""The Fortran proximal operators must agree with the NumPy reference.""" + +import numpy as np +import pytest +import reference + +import smopt + + +SHAPES = [(1, 1), (5, 1), (12, 3), (40, 7)] +GAMMAS = [0.0, 0.05, 0.5, 4.0] + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +@pytest.mark.parametrize("gamma", GAMMAS) +def test_prox_l1_matches_reference( + rng: np.random.Generator, n: int, p: int, gamma: float +) -> None: + x = rng.standard_normal((n, p)) + got = smopt.prox_l1(x, 0.7, gamma=gamma) + assert np.allclose(got, reference.prox_l1(x, 0.7, gamma=gamma)) + + +@pytest.mark.parametrize(("n", "p"), SHAPES) +@pytest.mark.parametrize("gamma", GAMMAS) +def test_prox_l21_matches_reference( + rng: np.random.Generator, n: int, p: int, gamma: float +) -> None: + x = rng.standard_normal((n, p)) + got = smopt.prox_l21(x, 0.7, gamma=gamma) + assert np.allclose(got, reference.prox_l21(x, 0.7, gamma=gamma)) + + +def test_prox_with_zero_weight_is_the_identity( + rng: np.random.Generator, +) -> None: + x = rng.standard_normal((10, 3)) + assert np.allclose(smopt.prox_l1(x, 1.0), x) + assert np.allclose(smopt.prox_l21(x, 1.0), x) + + +def test_prox_l1_soft_thresholds() -> None: + x = np.array([[-3.0, -0.5, 0.0, 0.5, 3.0]]) + got = smopt.prox_l1(x, 1.0, gamma=1.0) + assert np.allclose(got, [[-2.0, 0.0, 0.0, 0.0, 2.0]]) + + +def test_prox_l21_shrinks_whole_rows() -> None: + x = np.array([[3.0, 4.0], [0.3, 0.4]]) + got = smopt.prox_l21(x, 1.0, gamma=1.0) + # The first row has norm 5 and keeps its direction scaled by 4/5; + # the second has norm 0.5 and is annihilated. + assert np.allclose(got, [[2.4, 3.2], [0.0, 0.0]]) + + +def test_prox_l21_survives_a_zero_row() -> None: + x = np.zeros((3, 2)) + assert np.allclose(smopt.prox_l21(x, 1.0, gamma=1.0), 0.0) + + +@pytest.mark.parametrize("gamma", [0.1, 1.0]) +def test_prox_l1_solves_its_own_subproblem( + rng: np.random.Generator, gamma: float +) -> None: + """The prox must beat nearby points on its defining objective.""" + x = rng.standard_normal((8, 3)) + eta = 0.3 + y = smopt.prox_l1(x, eta, gamma=gamma) + + def value(z: np.ndarray) -> float: + return float( + np.sum((z - x) ** 2) / (2 * eta) + gamma * np.sum(np.abs(z)) + ) + + best = value(y) + for _ in range(20): + assert best <= value(y + 1e-3 * rng.standard_normal(y.shape)) + 1e-12