groupby-sum pads every group to the largest group size

[FBumann]

Note: This issue was drafted by an AI agent (Claude Code) during a profiling/investigation session with @FBumann. The measurements (memray C-level profiling, benchmarks, the 284× / 132 MB figures) and the git-history archaeology are real and reproducible, but the framing and proposed directions are a starting point for discussion — please sanity-check before acting on them.

LinearExpression.groupby(...).sum() pads every group to the largest group size along _term (via unstack(group_dim, fill_value=...) in expressions.py). For skewed groups — e.g. incidence/balance constraints where a few buses host most units — the result blows up: a 300k-element expression over zipf-distributed groups produces a ~12 GiB padded result that is >99% fill.

Two separable costs:

• Transient (the unstack/stack machinery making extra copies) — addressable in the current representation. Prototype: perf: scatter groupby-sum terms directly instead of unstacking FBumann/linopy#25 scatters terms directly into the result (~32× faster, ~2.6× less peak on the skewed case).
• Inherent (the padded result itself) — forced by the dense rectangular _term layout; only a long-format/sparse kernel removes it. See Umbrella: long-format / sparse _term kernel (dense-_term memory cluster) #756.

Real-world evidence: PyPSA already works around this

PyPSA splits the nodal-balance constraint into two separate constraints — strongly-meshed vs weakly-meshed buses — purely to contain this padding (pypsa/optimization/optimize.py, with its own comment "This reduces memory usage for large networks"):

meshed_buses = get_strongly_meshed_buses(n, threshold=45)   # buses with >45 attached components
define_nodal_balance_constraints(n, sns, buses=weakly_meshed_buses)             # small groups
define_nodal_balance_constraints(n, sns, buses=meshed_buses, suffix="-meshed")  # large groups

A single grouped balance over all buses would pad every bus to the largest hub's term count; bucketing by meshedness keeps each rectangle padded only to its own bucket max. This is the "eventually do a separation of short and long linear expressions" noted in the original groupby commit (PyPSA #557, 2023) — and it is actively maintained: PyPSA #1591 (2026) promoted meshed_thresholds to a tunable user parameter.

Note: on typical PyPSA networks the realistic group skew is small (max ≈ 8 generators/bus, ~2.7× padding), so groupby is not the build's peak allocation there — merge is (#749). The pathological skew matters for detailed unit-commitment / rooftop-PV aggregation and for the meshed hubs the split above exists to handle.

Related

• @/dot against a sparse matrix densifies the result to full _term #748 — @/dot against a sparse matrix densifies the same way (KVL).
• merge of ragged expressions is the peak allocation in PyPSA model builds #749 — merge of ragged expressions is the actual build peak.
• Umbrella: long-format / sparse _term kernel (dense-_term memory cluster) #756 — long-format/sparse _term kernel that subsumes all three.

Filing as a tracking note so this doesn't get lost.

[FBumann]

Note: Follow-up from the same profiling session (AI agent + @FBumann). Reproducible; framing is a starting point.

Env: pypsa 1.2.2.post1.dev2+g30e4ed0ef (git 30e4ed0e), linopy 0.7.0, xarray 2026.4.0, numpy 2.4.6, pandas 2.3.3, scipy 1.16.3, python 3.13.2, macOS-26.2 arm64; memray 1.19.3 (native traces); build-only (create_model, no solve).

1. The padding cost, measured by toggling PyPSA's own workaround. create_model(meshed_thresholds=[]) collapses the default [30,100,400,inf] bus-banding into a single nodal-balance constraint — i.e. turns the split OFF on a real build path. Synthetic radial net (no cycles → KVL trivial, isolates the balance): 120 buses, 19,900 generators skewed 8000 on one bus / 100 on each other, 24 snapshots → 480,456 vars.

config	build peak	total allocated
default banding (split on)	95 MB	0.70 GB
meshed_thresholds=[] (split off)	1.16 GB	2.41 GB

12× peak from the padding alone — the banding is load-bearing, not cosmetic.

2. Where the width is created — before merge, in groupby. Tracing _term: the final Bus-nodal_balance is name=120 × snapshot=24 × _term=8120 = 23.4M coeff cells. The merge that builds it (constraints.py:1139) decomposes as:

merge inputs (per component, after groupby().sum()):
   Generator block  _term = 8000   <- padded to the 8000-gen hub
   Line block       _term =  119
   Load block       _term =    1
merge output        _term = 8120   ( = 8000 + 119 + 1 )

8000 of 8120 terms exist before merge, created by groupby(bus).sum() padding every bus to the largest group (this issue). Merge along _term is faithful (out = Σ in; no extra padding) — it only materializes the already-padded blocks and broadcasts them across all 120 buses.

Fix location: upstream of merge. A sparse/long-format _term in groupby().sum() (or #740) shrinks the 8000-wide block at the source and merge inherits the saving; a merge-side change can't recover it, because the bloat is in merge's inputs.