Port the Stiefel manifold solvers to Fortran 77

.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

README.md

@@ -1,42 +1,115 @@
-# <Project Name>
+![smopt](https://raw.githubusercontent.com/eggzec/smopt/master/docs/assets/images/smopt-banner.png)
 
-**<Project Description>**
+# smopt
 
-[![Tests](https://github.com/eggzec/<Project Name>/actions/workflows/test.yml/badge.svg)](https://github.com/eggzec/<Project Name>/actions/workflows/test.yml)
-[![Documentation](https://github.com/eggzec/<Project Name>/actions/workflows/docs.yml/badge.svg)](https://github.com/eggzec/<Project Name>/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/<Project Name>/graph/badge.svg)](https://codecov.io/github/eggzec/<Project Name>)
-[![Quality Gate Status](https://sonarcloud.io/api/project_badges/measure?project=eggzec_<Project Name>&metric=alert_status)](https://sonarcloud.io/project/overview?id=eggzec_<Project Name>)
+[![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/<Project Name>.svg?label=PyPI%20downloads)](https://pypi.org/project/<Project Name>/)
-[![Python versions](https://img.shields.io/pypi/pyversions/<Project Name>.svg)](https://pypi.org/project/<Project Name>/)
+[![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.
+$$
 
-`<Project Description>`
+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 <Project Name>
+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 <Project Name>
+pip install smopt
 ```
 
-Requires Python 3.10+ and NumPy. No external runtime dependencies. See the
-[full installation guide](https://eggzec.github.io/<Project Name>/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/<Project Name>/theory/) — mathematical background, hierarchical basis, algorithms
-- [Quickstart](https://eggzec.github.io/<Project Name>/quickstart/) — runnable examples
-- [API Reference](https://eggzec.github.io/<Project Name>/api/) — class and function signature and arguments
-- [References](https://eggzec.github.io/<Project Name>/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
+(<https://stmopt.gitee.io/>). `smopt` re-implements its numerics in
+Fortran 77 behind the same solver interface.
 
 ## License
 

bin/build.py

@@ -78,7 +78,7 @@ def wheel():
 def clean():
     logger.debug("Starting cleanup ...")
 
-    run_command("uv pip uninstall <Project Name>")
+    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="<Project Name> Build Script")
+    parser = argparse.ArgumentParser(description="smopt Build Script")
     parser.add_argument(
         "mode",
         help="""Build mode:

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 <signature.pyf> <build-dir> [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())

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))
+```