Docs p1

README.md

@@ -1,72 +1,170 @@
 # within
 
-High-performance fixed-effects solver for econometric panel data. Solves
-**y = D x** where D is a sparse categorical design matrix, using
-iterative methods (preconditioned CG and right-preconditioned GMRES) with
-domain decomposition (Schwarz) preconditioners backed by approximate Cholesky
-local solvers.
+`within` provides high-performance solvers for projecting out high-dimensional fixed effects from regression problems.
 
-Designed for millions of observations with multiple high-dimensional fixed
-effects.
+By the Frisch-Waugh-Lovell theorem, estimating a regression of the form *y = Xβ + Dα + ε* reduces to a sequence of least-squares projections, one for y and one for each column of X, followed by a cheap regression fit on the resulting residuals. The projection step of solving the normal equations *D'Dx = D'z* is the computational bottleneck, which is the problem `within` is designed to solve.
 
-## Python quickstart
+`within`'s solvers are tailored to the structure of fixed effects problems, which can be represented as a graph (as first noted by Correia, 2016). Concretely, `within` uses iterative methods (preconditioned CG, right-preconditioned GMRES) with domain decomposition (Schwarz) preconditioners, backed by approximate Cholesky local solvers (Gao et al, 2025).
+
+## Installation
 
 Requires Python >= 3.9.
 
-```bash
-pip install within
-```
+## Python Quickstart
+
+`within`'s main user-facing function is `solve`. Provide a 2-D `uint32` array of category codes (one column per fixed-effect factor) and a response vector `y`. The solver finds x in the normal equations **D'D x = D'y**, where D is the sparse categorical design matrix.
 
 ```python
-from within import solve, CG, GMRES, MultiplicativeSchwarz
+from within import solve, solve_batch
 import numpy as np
 
-# Two factors with 500 levels each, 100k observations
-categories = [
-    np.random.randint(0, 500, size=100_000),
-    np.random.randint(0, 500, size=100_000),
-]
-y = np.random.randn(100_000)
-
-# Solve with Schwarz-preconditioned CG (default)
-result = solve(categories, y)
-print(f"converged={result.converged}  iters={result.iterations}  residual={result.residual:.2e}")
-
-# Solve with unpreconditioned CG
-result = solve(categories, y, CG(preconditioner=None))
-print(f"converged={result.converged}  iters={result.iterations}")
-
-# Solve with GMRES + multiplicative Schwarz
-result = solve(categories, y, GMRES(preconditioner=MultiplicativeSchwarz()))
-print(f"converged={result.converged}  iters={result.iterations}")
-
-# Weighted solve
-weights = np.random.exponential(1.0, size=100_000)
-result = solve(categories, y, weights=weights)
-print(f"converged={result.converged}  iters={result.iterations}")
+np.random.seed(1)
+n = 100_000
+fe = np.asfortranarray(np.column_stack([
+    np.random.randint(0, 500, n).astype(np.uint32),
+    np.random.randint(0, 200, n).astype(np.uint32),
+]))
+y = np.random.randn(n)
+
+result = solve(fe, y)                          # Schwarz-preconditioned CG
+result = solve(fe, y, weights=np.ones(n))      # weighted solve
+```
+
+### FWL regression example
+
+```python
+beta_true = np.array([1.0, -2.0, 0.5])
+X = np.random.randn(n, 3)
+y = X @ beta_true + np.random.randn(n)
+
+result = solve_batch(fe, np.column_stack([y, X]))
+y_tilde, X_tilde = result.demeaned[:, 0], result.demeaned[:, 1:]
+beta_hat = np.linalg.lstsq(X_tilde, y_tilde, rcond=None)[0]
+print(np.round(beta_hat, 4))  # [ 0.9982 -2.006   0.5005]
+```
+
+## Python API
+
+### High-level functions
+
+| Function | Description |
+|---|---|
+| `solve(categories, y, config?, weights?, preconditioner?)` | Solve a single right-hand side. Returns `SolveResult`. |
+| `solve_batch(categories, Y, config?, weights?, preconditioner?)` | Solve multiple RHS vectors in parallel. `Y` has shape `(n_obs, k)`. Returns `BatchSolveResult`. |
+
+`categories` is a 2-D `uint32` array of shape `(n_obs, n_factors)`. A `UserWarning` is emitted when a C-contiguous array is passed — use `np.asfortranarray(categories)` for best performance.
+
+### Persistent solver
+
+For repeated solves with the same design matrix, `Solver` builds the preconditioner once and reuses it.
+
+```python
+from within import Solver
+
+solver = Solver(fe)
+r = solver.solve(y)                            # reuses preconditioner
+r = solver.solve_batch(np.column_stack([y, X]))
+
+precond = solver.preconditioner()              # picklable
+solver2 = Solver(fe, preconditioner=precond)   # skip re-factorization
 ```
 
-### Solver methods
+| Property / Method | Description |
+|---|---|
+| `Solver(categories, config?, weights?, preconditioner?)` | Build solver. Factorizes the preconditioner at construction. |
+| `.solve(y)` | Solve a single RHS. Returns `SolveResult`. |
+| `.solve_batch(Y)` | Solve multiple RHS columns in parallel. Returns `BatchSolveResult`. |
+| `.preconditioner()` | Return the built `FePreconditioner` (picklable), or `None`. |
+
+
+### Solver configuration
 
 | Class | Description |
 |---|---|
-| `CG(tol, maxiter, preconditioner, operator)` | Conjugate gradient on the normal equations. Default is additive Schwarz with the implicit operator. |
-| `GMRES(tol, maxiter, restart, preconditioner, operator)` | Right-preconditioned GMRES on the normal equations. |
+| `CG(tol=1e-8, maxiter=1000, operator=Implicit)` | Conjugate gradient on the normal equations. Default solver. |
+| `GMRES(tol=1e-8, maxiter=1000, restart=30, operator=Implicit)` | Right-preconditioned GMRES. |
 
 ### Preconditioners
 
 | Class | Description |
 |---|---|
-| `AdditiveSchwarz(...)` | Additive one-level Schwarz. Use with `CG` or `GMRES`. |
-| `MultiplicativeSchwarz(...)` | Multiplicative one-level Schwarz. Use with `GMRES` only. |
+| `AdditiveSchwarz(local_solver?)` | Additive one-level Schwarz. Works with CG and GMRES. |
+| `MultiplicativeSchwarz(local_solver?)` | Multiplicative one-level Schwarz. GMRES only. |
+| `Preconditioner.Off` | Disable preconditioning. |
 
-## Rust quickstart
+Pass `None` (the default) to use additive Schwarz with the default local solver.
 
-```bash
-cargo add within
+### Local solver configuration (advanced)
+
+| Class | Description |
+|---|---|
+| `SchurComplement(approx_chol?, approx_schur?, dense_threshold=24)` | Schur complement reduction with approximate Cholesky on the reduced system. Default local solver. |
+| `FullSddm(approx_chol?)` | Full bipartite SDDM factorized via approximate Cholesky. |
+| `ApproxCholConfig(seed=0, split=1)` | Approximate Cholesky parameters. |
+| `ApproxSchurConfig(seed=0, split=1)` | Approximate Schur complement sampling parameters. |
+
+### Result types
+
+**`SolveResult`**: `x` (coefficients), `demeaned` (residuals), `converged`, `iterations`, `residual`, `time_total`, `time_setup`, `time_solve`.
+
+**`BatchSolveResult`**: Same fields, with `converged`, `iterations`, `residual`, and `time_solve` as lists (one entry per RHS).
+
+## Rust API
+
+```rust
+use ndarray::Array2;
+use within::{solve, SolverParams, KrylovMethod, Preconditioner, LocalSolverConfig};
+
+let categories = /* Array2<u32> of shape (n_obs, n_factors) */;
+let y: &[f64] = /* response vector */;
+
+// Default: CG + additive Schwarz
+let r = solve(categories.view(), &y, None, &SolverParams::default(), None)?;
+assert!(r.converged);
+
+// GMRES + multiplicative Schwarz
+let params = SolverParams {
+    krylov: KrylovMethod::Gmres { restart: 30 },
+    ..SolverParams::default()
+};
+let precond = Preconditioner::Multiplicative(LocalSolverConfig::default());
+let r = solve(categories.view(), &y, None, &params, Some(&precond))?;
 ```
 
-See [`crates/within/README.md`](crates/within/README.md) for Rust API usage and architecture details.
+Persistent solver — build once, solve many:
+
+```rust
+use within::Solver;
+
+let solver = Solver::new(categories.view(), None, &SolverParams::default(), None)?;
+let r1 = solver.solve(&y)?;
+let r2 = solver.solve(&another_y)?;  // reuses preconditioner
+```
+
+| Type | Variants / Fields |
+|---|---|
+| `SolverParams` | `krylov: KrylovMethod`, `operator: OperatorRepr`, `tol: f64`, `maxiter: usize` |
+| `KrylovMethod` | `Cg` (default), `Gmres { restart }` |
+| `OperatorRepr` | `Implicit` (default, matrix-free D'WD), `Explicit` (CSR Gramian) |
+| `Preconditioner` | `Additive(LocalSolverConfig)`, `Multiplicative(LocalSolverConfig)` |
+| `LocalSolverConfig` | `SchurComplement { approx_chol, approx_schur, dense_threshold }`, `FullSddm { approx_chol }` |
+
+### Lower-level types
+
+The crate exposes its internals for advanced use:
+
+| Module | Key types |
+|---|---|
+| `observation` | `FactorMajorStore`, `ArrayStore`, `ObservationStore` trait |
+| `domain` | `WeightedDesign`, `FixedEffectsDesign`, `Subdomain` |
+| `operator` | `Gramian` (CSR), `GramianOperator` (implicit), `DesignOperator`, `build_schwarz`, `FeSchwarz` |
+| `solver` | `Solver<S: ObservationStore>` — generic persistent solver |
+
+### Feature flags
+
+| Feature | Default | Effect |
+|---|---|---|
+| `ndarray` | yes | Enables `from_array` constructors for `ndarray::ArrayView2` interop. |
 
 ## Project structure
 
@@ -76,7 +174,7 @@ crates/
   within/            Core fixed-effects solver (observation stores, domains, operators, orchestration)
   within-py/         PyO3 bridge (cdylib → within._within)
 python/within/       Python package re-exporting the Rust extension
-benchmarks/          Python benchmark framework (18 suites)
+benchmarks/          Python benchmark framework
 ```
 
 ## Development
@@ -99,18 +197,7 @@ MIT
 
 ## References
 
-- Correia, S. (2017). Linear Models with High-Dimensional Fixed Effects: An
-  Efficient and Feasible Estimator.
-  https://hdl.handle.net/10.2139/ssrn.3129010
-
-- Gaure, S. (2013). OLS with Multiple High Dimensional Category Variables.
-  *Computational Statistics & Data Analysis*, 66, 2--17.
-
-- Gao, Y., Kyng, R. & Spielman, D. A. (2025). AC(k): Robust Solution of
-  Laplacian Equations by Randomized Approximate Cholesky Factorization.
-  *SIAM Journal on Scientific Computing*.
-
-- Xu, J. (1992). Iterative Methods by Space Decomposition and Subspace
-  Correction. *SIAM Review*, 34(4), 581--613.
-
+- Correia, Sergio. "A feasible estimator for linear models with multi-way fixed effects." *Preprint* at http://scorreia.com/research/hdfe.pdf (2016).
+- Gao, Y., Kyng, R. & Spielman, D. A. (2025). AC(k): Robust Solution of Laplacian Equations by Randomized Approximate Cholesky Factorization. *SIAM Journal on Scientific Computing*.
 - Toselli & Widlund (2005). *Domain Decomposition Methods — Algorithms and Theory*. Springer.
+- Xu, J. (1992). Iterative Methods by Space Decomposition and Subspace Correction. *SIAM Review*, 34(4), 581--613.

crates/within/README.md

@@ -16,66 +16,60 @@ cargo add within
 ## Quick example
 
 ```rust
-use within::{
-    solve, GmresPrecond, LocalSolverConfig, OperatorRepr, SolverMethod, SolverParams,
-};
+use ndarray::Array2;
+use within::{solve, SolverParams, KrylovMethod, Preconditioner, LocalSolverConfig};
 
 // Two factors: 100 levels each, 10 000 observations
-let factor_0: Vec<u32> = (0..10_000).map(|i| (i % 100) as u32).collect();
-let factor_1: Vec<u32> = (0..10_000).map(|i| (i / 100) as u32).collect();
-let y: Vec<f64> = (0..10_000).map(|i| i as f64 * 0.01).collect();
+let n_obs = 10_000usize;
+let mut categories = Array2::<u32>::zeros((n_obs, 2));
+for i in 0..n_obs {
+    categories[[i, 0]] = (i % 100) as u32;
+    categories[[i, 1]] = (i / 100) as u32;
+}
+let y: Vec<f64> = (0..n_obs).map(|i| i as f64 * 0.01).collect();
 
 // Solve with the default solver: CG + additive Schwarz + implicit operator
-let result = solve(
-    &[factor_0.clone(), factor_1.clone()],
-    &[100, 100],
-    &y,
-    None,
-    &SolverParams::default(),
-).expect("solve should succeed");
+let result = solve(categories.view(), &y, None, &SolverParams::default(), None)
+    .expect("solve should succeed");
 assert!(result.converged);
-println!("default CG converged in {} iterations", result.iterations);
+println!("CG converged in {} iterations", result.iterations);
 
-// Solve with GMRES + multiplicative Schwarz on the implicit operator
+// GMRES with multiplicative Schwarz
 let params = SolverParams {
-    method: SolverMethod::Gmres {
-        preconditioner: Some(GmresPrecond::Multiplicative(LocalSolverConfig::gmres_default())),
-        operator: OperatorRepr::Implicit,
-        restart: 30,
-    },
-    tol: 1e-8,
-    maxiter: 1000,
+    krylov: KrylovMethod::Gmres { restart: 30 },
+    ..SolverParams::default()
 };
-let result = solve(
-    &[factor_0, factor_1],
-    &[100, 100],
-    &y,
-    None,
-    &params,
-).expect("solve should succeed");
+let precond = Preconditioner::Multiplicative(LocalSolverConfig::default());
+let result = solve(categories.view(), &y, None, &params, Some(&precond))
+    .expect("solve should succeed");
 assert!(result.converged);
 println!("GMRES converged in {} iterations", result.iterations);
 ```
 
+## Feature flags
+
+| Feature   | Default | Effect                                                                 |
+|-----------|---------|------------------------------------------------------------------------|
+| `ndarray` | yes     | Enables `from_array` constructors on observation stores for interop with `ndarray::ArrayView2`. |
+
 ## Architecture
 
 The crate is organized in four layers:
 
-1. **`observation`** -- Per-observation factor levels and weights via
+1. **`observation`** — Per-observation factor levels and weights via
    `FactorMajorStore` and the `ObservationStore` trait.
 
-2. **`domain`** -- Domain decomposition. `WeightedDesign` wraps a store with
+2. **`domain`** — Domain decomposition. `WeightedDesign` wraps a store with
    factor metadata; `build_local_domains` constructs factor-pair subdomains
    with partition-of-unity weights for the Schwarz preconditioner.
 
-3. **`operator`** -- Linear algebra primitives. `Gramian` (explicit CSR) and
+3. **`operator`** — Linear algebra primitives. `Gramian` (explicit CSR) and
    `GramianOperator` (implicit D^T W D) for the normal equations; Schwarz
    preconditioner builders that wire approximate Cholesky local solvers into
    the generic `schwarz-precond` framework.
 
-4. **`orchestrate`** -- End-to-end solve entry points (`solve`,
-   `solve_normal_equations`) with
-   typed configuration (`SolverParams`, `SolverMethod`, `GmresPrecond`,
+4. **`orchestrate`** — End-to-end solve entry points (`solve`, `solve_batch`)
+   with typed configuration (`SolverParams`, `KrylovMethod`, `Preconditioner`,
    `OperatorRepr`).
 
 ## License
@@ -84,15 +78,7 @@ MIT
 
 ## References
 
-- Correia, S. (2017). Linear Models with High-Dimensional Fixed Effects: An
-  Efficient and Feasible Estimator. *Working paper*.
-  https://hdl.handle.net/10.2139/ssrn.3129010
-
-- Gaure, S. (2013). OLS with Multiple High Dimensional Category Variables.
-  *Computational Statistics & Data Analysis*, 66, 2--17.
-  https://doi.org/10.1016/j.csda.2013.03.024
-
-- Spielman, D. A. & Teng, S.-H. (2014). Nearly Linear Time Algorithms for
-  Preconditioning and Solving Symmetric, Diagonally Dominant Linear Systems.
-  *SIAM Journal on Matrix Analysis and Applications*, 35(3), 835--885.
-  https://doi.org/10.1137/090771430
+- Correia, Sergio. "A feasible estimator for linear models with multi-way fixed effects." *Preprint* at http://scorreia.com/research/hdfe.pdf (2016).
+- Gao, Y., Kyng, R. & Spielman, D. A. (2025). AC(k): Robust Solution of Laplacian Equations by Randomized Approximate Cholesky Factorization. *SIAM Journal on Scientific Computing*.
+- Toselli & Widlund (2005). *Domain Decomposition Methods — Algorithms and Theory*. Springer.
+- Xu, J. (1992). Iterative Methods by Space Decomposition and Subspace Correction. *SIAM Review*, 34(4), 581--613.