Benchmark comparison vs `polar-high` (Polars-backed LP/MIP eDSL on HiGHS) — thoughts?

[MaykThewessen]

A new modelling library, polar-high, was released recently by Juha Kiviluoma (Nodal-Tools Oy). It targets the same niche as linopy: an indexed LP/MIP eDSL on top of HiGHS, but uses a Polars-backed kernel instead of xarray.

The docs include a reproducible build / solve / memory comparison between polar-high, linopy, and Pyomo on indexed linear programs, all using HiGHS as the solver. The dense-LP test follows linopy's own benchmark, restricted to the linear case so polar-high can run it directly.

Dense LP — build-only headline (single-thread, HiGHS time-limited to ~1 µs to isolate modelling-layer cost)

At N=3000 (~18M variables):

Tool	Total time	Peak memory
polar-high (regular)	87 s	16.4 GB
polar-high (save_memory)	131 s	15.1 GB
linopy (io_api="lp")	300 s	23.7 GB

Network LP — irregular topology (N nodes, 5N edges, T=168, build-only, 1 thread)

At N=10000 (~8.4M variables):

Tool	Total time	Peak memory
polar-high (regular)	51 s	11.5 GB
polar-high (save_memory)	94 s	12.2 GB
linopy (io_api="lp")	165 s	14.1 GB

Replication: linopy's original benchmark format (3-panel, full HiGHS solve included)

The replication is faithful: same dense LP, same closed-form optimum, HiGHS solves the same sparse matrix. linopy's published numbers reproduce within run-to-run noise.

Summary

Reported ~2.3-3.5x faster than linopy on the build path with ~30% lower peak memory. The Polars-backed lazy/columnar internals are the architectural source of the gap; linopy is xarray/dense-numpy by default.

Questions for the maintainers and community

1. Methodology — The linopy column uses io_api="lp" because the benchmark's time_limit short-circuit interacts cleanly with that path. Would io_api="direct" materially change the picture on the build column?
2. Architectural direction — Issues Polars support #510 (Polars support) and lp-polars IO API writes different models to lp which are incorrect #484 (lp-polars IO bugs) both touched a Polars path. Given the gap shown here, is a Polars-backed write path inside linopy worth revisiting? Or is the architectural distance from xarray too large?
3. Cross-tool collaboration — Anyone here in contact with Juha / Nodal-Tools? A shared LP/MIP benchmark suite across linopy / polar-high / Pyomo (and ideally JuMP) would be useful for PyPSA-Eur-class workloads.

Not affiliated with Nodal-Tools, just a downstream PyPSA-Eur user who finds the numbers interesting.

[FBumann]

@MaykThewessen Thanks for raising this here. I was not aware of this project. I think polars is quite promising as a LP/MILP backend.

[MaykThewessen]

Thanks, i saw it in my Linkedin feed on: linkedin.com/messaging/thread/2-NmZkMGI5NjYtOTQ2MS00YTlkLWFkMzMtYzEyN2FhYWRkMjk0XzEwMA==/

[FabianHofmann]

yes I think so too. definitely want to have a look at their project. @MaykThewessen yes io_api="direct" will make a difference. if you ran these, could you quickly check out the impact?

[MaykThewessen]

Follow-up: tested the in-tree CSRConstraint / freeze_constraints path (#630)

Worth noting linopy already has an in-tree memory-lean path from #630 (CSRConstraint, the freeze_constraints flag), which is arguably the more relevant comparison point than a full Polars rewrite. I tested it on a real PyPSA DC-OPF network (NL + neighbouring bidding zones, 168 h slice, ~1.78M constraint rows) on linopy 0.7.0.post1.dev19+g37af4ba12.

Plumbing through PyPSA works out of the box:

n.optimize(model_kwargs={"freeze_constraints": True})

PyPSA's create_model forwards **kwargs to Model(**kwargs), and the constraint backing flips from Constraint to CSRConstraint as expected.

Measured on create_model (one variant per process for a clean peak):

	constraint storage	process peak RSS	build
default (xarray)	242.5 MB	4401 MB	2.19 s
freeze_constraints=True (CSR)	165.3 MB	4376 MB	1.93 s

So on this real (sparse, non-padded) network LP the CSR path is ~32% leaner on constraint storage, not the ~90% headline from the padded dense benchmarks. Process peak RSS barely moved (~0.6%) here because constraint storage is only ~5% of the process footprint in my proxy (which loaded an already-solved network, inflating the baseline; I did not isolate the pure pre-solve build peak).

Two takeaways for the questions above:

1. Re Polars backing (Q2): there is already a CSR-backed lean path in linopy. Before concluding a Polars rewrite is needed, it would be worth re-running the polar-high comparison with io_api="direct" and freeze_constraints=True together, since that is the closest apples-to-apples lean configuration linopy currently offers.
2. The realistic memory win from CSR on sparse network LPs looks closer to ~30% than the dense-test headline, which narrows (but does not close) the gap to polar-high's reported ~30% lower peak.

Happy to share the small measurement script if useful.

[FBumann]

Note

AI-written (Claude Code, prompted by @FBumann).

Another concrete data point for the dense-vs-sparse kernel question, surfaced while implementing #753 (multi-key groupby([names])).

LinearExpression.groupby([names]).sum() returns one dimension per key — which in dense xarray is a full cartesian grid, so every absent key combination becomes a real fill cell. For a sparse/correlated crossing this explodes. Measured with memray, a diagonal crossing (N=1000 observed combos):

grouping	output	peak memory
groupby([names])	{period:1000, season:1000} dense grid	33.3 MB
groupby(df)	{group:1000} MultiIndex, observed-only	0.33 MB

~100×, scaling as N. The difference is entirely the final densification.

Takeaway for the kernel decision: "separate dims + compact" is impossible in dense xarray — dims are dense, so a grid is the only separate-dims representation. Staying compact requires either a stacked MultiIndex (the DataFrame-grouper path) or a genuinely sparse backing: sparse.COO duck-arrays under xarray, or a long-format/polars kernel. Under a sparse/long kernel this blow-up, the matmul zero-padding (#748), and the merge padding (#749) all dissolve together — none would densify.

So groupby is one more op where the dense _term/dense-dims model forces materialization a sparse kernel would avoid. Evidence on the "dense xarray ceiling" side.