You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Fast and accurate weighted least squares meshless interpolator for Python.
17
-
18
16
We use [semantic versioning](https://semver.org/).
19
17
20
18
For my stance on AI contributions, see the [collaboration guidelines](https://github.com/Technologicat/substrate-independent/blob/main/collaboration.md).
@@ -24,64 +22,88 @@ For my stance on AI contributions, see the [collaboration guidelines](https://gi
24
22
25
23
## Introduction
26
24
27
-
WLSQM (Weighted Least SQuares Meshless) constructs a piecewise polynomial surrogate model on scattered data: given scalar values at a point cloud in 1D, 2D, or 3D, it fits a local polynomial (up to 4th order) in the neighborhood of each query point, by weighted least squares. From the surrogate you can read off the function value and any derivative up to the polynomial order, at the fit origin or anywhere inside the local neighborhood.
25
+
The `wlsqm` Python library is a fast and accurate weighted least squares meshless interpolator and differentiator.
26
+
27
+
The WLSQM (Weighted Least SQuares Meshless) method constructs a piecewise polynomial surrogate model on scattered data. Given scalar values at a point cloud in 1D, 2D, or 3D, this library fits a local polynomial (up to 4th order) in the neighborhood of each data point, by weighted least squares. From the surrogate model, you can read off the function value and any derivative up to the polynomial order, at the fit origin (each data point) or anywhere inside the local neighborhood.
28
28
29
29
Use cases:
30
30
31
31
-**Numerical differentiation** of data known only as values at discrete points. Applies to explicit algorithms for IBVPs: compute a gradient or Laplacian on a scattered point cloud without needing a mesh.
32
32
-**Response-surface modeling** on an unstructured set of design points.
33
-
-**Smoothing noisy function values** — the averaging effect of the least-squares fit denoises first derivatives significantly. (Second derivatives are more sensitive to noise; for those, the robust recipe is to run WLSQM once to recover first derivatives on the neighborhood, then run it again on those to get second derivatives.)
33
+
-**Smoothing noisy function values** — the averaging effect of the least-squares fit denoises first derivatives significantly.
34
+
- Second derivatives are more sensitive to noise; for those, the robust recipe is to run WLSQM once to recover first derivatives on the neighborhood, then run it again on those to get second derivatives.
34
35
35
36
No grid or mesh is needed. The only restriction on geometry is non-degeneracy — e.g. 2D points must not all fall on the same 1D line.
36
37
37
38
### Not a Taylor series
38
39
39
-
Despite the storage layout looking like a Taylor expansion (slots for `f`, `df/dx`, `df/dy`, `d²f/dx²`, ..., indexed by monomial degree), WLSQM is **not** evaluating a Taylor series. The polynomial is fitted over a local neighborhood of data points by weighted least squares; the coefficients come from a linear solve, not from differentiating an analytic function at the origin. The error behavior is correspondingly much better than Taylor truncation error would predict, thanks to the averaging effect of the least-squares fit.
40
+
Despite the storage layout looking like the coefficients of a Taylor expansion (slots for `f`, `df/dx`, `df/dy`, `d²f/dx²`, ..., indexed by monomial degree), the WLSQM method is **not** derived from a Taylor series.
41
+
42
+
The polynomial is fitted over a local neighborhood of data points by weighted least squares. The coefficients come from a linear solve, not from differentiating an analytic function at the origin. The error behavior is correspondingly much better than Taylor truncation error would predict, thanks to the averaging effect of the least-squares fit.
40
43
41
-
This method appears in the literature under several names — MLS (moving least squares), WLSQM (weighted least squares meshless), "diffuse approximation." These are essentially the same idea with varying weighting function choices.
44
+
In a sense, we may even say that there *is no underlying function* that is being modeled. The only thing that exists before the fitting are the data values on a set of points. Whether or not these data values come from a differentiable function — or are just arbitrary numerical data — is immaterial.
42
45
43
-
### Academic reference
46
+
This method appears in the literature under several names: MLS (moving least squares), WLSQM (weighted least squares meshless/meshfree), "diffuse approximation." These are essentially the same idea with varying weighting function choices.
47
+
48
+
### Origin of this variant of the method
44
49
45
50
This is an independent implementation of the algorithm described (in the 2nd-order 2D case) in section 2.2.1 of Hong Wang (2012), *Evolutionary Design Optimization with Nash Games and Hybridized Mesh/Meshless Methods in Computational Fluid Dynamics*, Jyväskylä Studies in Computing 162, University of Jyväskylä. [ISBN 978-951-39-5007-1 (PDF)](http://urn.fi/URN:ISBN:978-951-39-5007-1)
46
51
47
-
Full theory for the generalized version (including the case of unknown function values, and the accuracy analysis): `doc/wlsqm_gen.pdf` in this repository.
52
+
Full derivation of the generalized version (including the case of unknown function values, and the accuracy analysis) can be found in this repository, in [`doc/wlsqm_gen.pdf`](doc/wlsqm_gen.pdf).
48
53
49
54
50
55
## Features
51
56
52
57
-**1D / 2D / 3D point clouds**, polynomials of order 0 through 4.
53
-
-**Any derivative** up to the polynomial order is available: at the fit origin as a DOF of the linear solve, and at any point inside the neighborhood via the interpolator. Differentiation of the basis polynomials is hardcoded for speed.
54
-
-**Knowns.** At the fit origin, any subset of `{f, ∂f/∂x, ∂f/∂y, ∂²f/∂x², …}` can be marked as known. The known DOFs are eliminated algebraically, shrinking the linear system to just the unknowns. The function value itself may be unknown — useful e.g. for Neumann boundary conditions in a PDE solver.
55
-
-**Weighting methods.**`WEIGHT_UNIFORM` gives the best overall fit of function values. `WEIGHT_CENTER` emphasizes points near the fit origin, improving derivative estimates there at the cost of fit quality far from the origin.
58
+
-**Any derivative** up to the polynomial order is available: at the fit origin as a DOF of the linear solve, and at any point inside the neighborhood via the interpolator.
59
+
- Differentiation of the basis polynomials is hardcoded for speed.
60
+
-**Knowns.** At the fit origin, any subset of `{f, ∂f/∂x, ∂f/∂y, ∂²f/∂x², …}` can be marked as known.
61
+
- The known DOFs are eliminated algebraically, shrinking the linear system to just the unknowns.
62
+
- The function value itself may be unknown at some of the data points — this is useful e.g. for Neumann boundary conditions in a PDE solver.
63
+
-**Weighting methods.**
64
+
-`WEIGHT_UNIFORM` weights all data points equally in each local neighborhood.
65
+
-`WEIGHT_CENTER` emphasizes points near the fit origin (each data point), improving derivative estimates there at the cost of fit quality far from the origin (of that neighborhood).
56
66
-**Sensitivity data.** Solution DOFs can be optionally differentiated w.r.t. the input function values, so you know how the output moves when the input wiggles.
57
-
-**Expert mode with separate prepare and solve stages.** If many fits share the same geometry but have different function-value data, `ExpertSolver.prepare()` generates and LU-factors the problem matrices once, and `solve()` can then be called many times with different `fk`. This is the fast path for time-stepping an IBVP over a fixed point cloud.
58
-
-**Parallel across independent local problems.**`fit_*_many_parallel` and `ExpertSolver(ntasks=N)` run independent fits across OpenMP threads. The linear-solver step is also parallel.
67
+
-**Expert mode with separate prepare and solve stages.**
68
+
- If many fits share the same geometry but have different function-value data, `ExpertSolver.prepare()` generates and LU-factors the problem matrices once, and `solve()` can then be called many times with different `fk`.
69
+
- This is the fast path for time-stepping an IBVP over a fixed point cloud.
70
+
-**Parallel across independent local problems.**
71
+
-`fit_*_many_parallel` and `ExpertSolver(ntasks=N)` run independent fits across OpenMP threads. The linear-solver step is also parallel across the problem instances.
59
72
-**Speed.**
60
73
- Performance-critical code is in Cython with the GIL released.
61
74
- LAPACK is called directly via [SciPy's Cython-level bindings](https://docs.scipy.org/doc/scipy/reference/linalg.cython_lapack.html); no GIL round-trip for the solver loop.
62
-
- The polynomial evaluator uses a symmetric Horner-like form with fused multiply-add (`fma`). See [`wlsqm/fitter/polyeval.pyx`](wlsqm/fitter/polyeval.pyx).
63
75
-**Accuracy.**
64
-
- Problem matrices are preconditioned by a symmetry-preserving iterative scaling (Ruiz 2001) before LU factorization, which is critical for high-order fits.
76
+
- Problem matrices are preconditioned by a symmetry-preserving iterative scaling (Ruiz, 2001) before LU factorization, which is critical for high-order fits.
77
+
- Reference: Daniel Ruiz. 2001. *A Scaling Algorithm to Equilibrate Both Rows and Columns Norms in Matrices*. Report RAL-TR-2001-034.
78
+
- The polynomial evaluator uses a custom symmetric form that is a natural multidimensional extension of the standard Horner form.
79
+
- Evaluation uses fused multiply-add (`fma`), thus rounding only once per accumulation step.
80
+
- For details, see [`wlsqm/fitter/polyeval.pyx`](wlsqm/fitter/polyeval.pyx).
65
81
- Optional iterative refinement inside the solver mitigates roundoff further.
66
-
- FMA in the polynomial evaluator rounds only once per accumulation step.
67
82
68
83
69
84
## Examples
70
85
71
-
Minimal 2D fit:
86
+
Minimal example using a manufactured solution to produce the input data: fit `f(x,y) = 1 + 2x + 3y + 4xy + 5x² + 6y²` on a scattered point cloud centered at the origin.
72
87
73
88
```python
74
89
import numpy as np
75
90
import wlsqm
76
91
77
-
# Fit f(x,y) = 1 + 2x + 3y + 4xy + 5x² + 6y² on a scattered point cloud
78
-
# centered at the origin.
92
+
# seed the RNG for reproducibility
79
93
rng = np.random.default_rng(42)
94
+
95
+
# point cloud
80
96
xk = rng.uniform(-1.0, 1.0, size=(30, 2))
97
+
98
+
# data values on the point cloud
81
99
fk =1+2*xk[:,0] +3*xk[:,1] +4*xk[:,0]*xk[:,1] +5*xk[:,0]**2+6*xk[:,1]**2
82
100
101
+
# point(s) where to evaluate the fit [[x0, y0], [x1, y1], ...]
83
102
xi = np.array([0.0, 0.0])
84
-
fi = np.zeros(wlsqm.number_of_dofs(2, 2))
103
+
104
+
# output array for the fit
105
+
fi = np.zeros(wlsqm.number_of_dofs(2, 2)) # (number_of_dimensions, polynomial_order)
106
+
85
107
wlsqm.fit_2D(
86
108
xk=xk, fk=fk, xi=xi, fi=fi,
87
109
sens=None, do_sens=False,
@@ -90,13 +112,16 @@ wlsqm.fit_2D(
90
112
debug=False,
91
113
)
92
114
93
-
# fi now holds the partial derivatives at (0, 0) in the order
94
-
# F, X, Y, X2, XY, Y2. For this exact polynomial, the fit recovers
95
-
# [1, 2, 3, 10, 4, 12] to machine precision.
115
+
# fi now holds the partial derivatives at (0, 0).
116
+
# The ordering is [f, ∂f/∂x, ∂f/∂y, ∂²f/∂x², ∂²f/∂x∂y, ∂²f/∂y²].
96
117
print(fi)
97
118
```
98
119
99
-
For a comprehensive tour of the API, see [`examples/wlsqm_example.py`](examples/wlsqm_example.py). For a minimal `ExpertSolver` example that demonstrates the prepare/solve separation, see [`examples/expertsolver_example.py`](examples/expertsolver_example.py).
120
+
For this exact polynomial, the fit recovers `[1, 2, 3, 10, 4, 12]` to machine precision.
121
+
122
+
For a comprehensive tour of the API, see [`examples/wlsqm_example.py`](examples/wlsqm_example.py).
123
+
124
+
For a minimal `ExpertSolver` example that demonstrates the prepare/solve separation, see [`examples/expertsolver_example.py`](examples/expertsolver_example.py).
100
125
101
126
102
127
## Installation
@@ -150,7 +175,7 @@ cd python-wlsqm
150
175
pdm config venv.in_project true
151
176
pdm use 3.14 # or whichever Python you prefer
152
177
pdm install # creates .venv, installs dev deps
153
-
export PATH="$(pwd)/.venv/bin:$PATH"# meson and ninja must be on PATH
178
+
export PATH="$(pwd)/.venv/bin:$PATH"#or $(pdm venv activate) — meson and ninja must be on PATH
-**API:** docstrings in [`wlsqm/fitter/simple.pyx`](wlsqm/fitter/simple.pyx) (simple API) and [`wlsqm/fitter/expert.pyx`](wlsqm/fitter/expert.pyx) (`ExpertSolver`).
173
207
-**Examples:**[`examples/wlsqm_example.py`](examples/wlsqm_example.py) for the full tour, [`examples/expertsolver_example.py`](examples/expertsolver_example.py) for `ExpertSolver` specifically.
174
208
-**Theory PDFs** in [`doc/`](doc/):
175
-
-[`wlsqm_gen.pdf`](doc/wlsqm_gen.pdf) — the generalized version (including the case of unknown function values), the accuracy analysis, and why WLSQM works.
176
-
-[`wlsqm.pdf`](doc/wlsqm.pdf) — older writeup for the pure-Python version originally written for FREYA, plus the sensitivity calculation.
209
+
-[`wlsqm_gen.pdf`](doc/wlsqm_gen.pdf) — the generalized version implemented in `wlsqm`, including the case of unknown function values, the accuracy analysis, and why WLSQM works.
210
+
-[`wlsqm.pdf`](doc/wlsqm.pdf) — older writeup, plus the sensitivity calculation.
177
211
-[`eulerflow.pdf`](doc/eulerflow.pdf) — application example in compressible flow, with a cleaner presentation of the original version.
178
212
179
213
See [TODO.md](TODO.md) for known gaps in the theory PDFs.
0 commit comments