StruMPL: Multi-task Dense Regression under Disjoint Partial Supervision and MNAR Labels
arXiv:2605.19931v1 [cs.CV] 19 May 2026
Reza M. Asiyabi∗ School of Geosciences, University of Edinburgh, UK National Centre for Earth Observation (NCEO), UK [email protected] Juan Alberto Molina-Valero Department of Spatial Sciences, Faculty of Environmental Sciences Czech University of Life Sciences Prague, Praha, Czech Republic [email protected] The SEOSAW Partnership School of Geosciences, University of Edinburgh, UK Steven Hancock School of Geosciences, University of Edinburgh, UK National Centre for Earth Observation (NCEO), UK [email protected] Casey M. Ryan School of Geosciences, University of Edinburgh, UK [email protected]
Abstract Estimating forest aboveground biomass (AGB) from Earth observation combines two structurally incompatible label sources: spaceborne lidar provides canopy structure at millions of locations but no biomass estimate, and ground-based plots provide biomass at thousands of biased locations but no metrics of structure. No single training sample carries labels for all target variables, plot labels are missing not at random (MNAR), and biomass is linked to the structural variables by known but biome-specific allometric laws. We formalise this as multi-task dense regression under heterogeneous disjoint partial supervision with MNAR labels and inter-task physical constraints, and propose StruMPL to address it jointly. A shared encoder feeds per-variable regression, imputation, and propensity heads for spatial MNAR correction, and a learnable physics module that evaluates the inter-task constraint on the model’s own predictions at every pixel. The supervised loss uses an Augmented IPW (AIPW) pseudo-outcome with stop-gradients on the propensity and on the imputation baseline; we show analytically and empirically that both are necessary for joint optimisation to recover IPW-weighted stationary points while keeping the loss bounded. On two ecologically distinct biomes, StruMPL outperforms ablation variants and the closest published method on AGB RMSE and bias, with a stratified analysis showing AIPW reduces high-AGB bias by ∼54%. ∗ https://reza-asiyabi.github.io/
Preprint.
1
Introduction
Aboveground biomass (AGB) is a major carbon stock in forests [Requena Suarez et al., 2019] and a required input for national greenhouse gas inventories and carbon offset markets [IPCC Task Force on National Greenhouse Gas Inventories (TFI), 2019]. Measuring AGB directly is destructive and infeasible at scale, so field ecologists measure structural attributes (e.g. height, cover and tree diameter) and estimate AGB through regionally calibrated allometric equations [Chave et al., 2014, Feldpausch et al., 2012]. Earth observation (EO) makes it possible to extend this indirect measurement chain to entire continents: spaceborne lidar (NASA’s GEDI [Dubayah et al., 2020]) provides canopy height and cover at millions of locations, and ground-based forest inventories provide AGB along with stem density and (sometimes) wood density [SEOSAW-partnership, 2021, Ministerio para la Transición Ecológica y el Reto Demográfico, 2024]. A model that jointly predicts all forest variables from the same input is desirable as: firstly, it would allow the structural and biomass predictions share a representation, and secondly it can produce maps that are bio-physically consistent with the allometric law that connects them. The two label sources, however, are structurally incompatible. GEDI footprints carry height and cover but no AGB; field plots carry AGB but no canopy structure at GEDI’s scale. No single training sample has labels for all target variables, and the source providing the most important target (AGB) accounts for less than 1% of the training data. Worse, field plots are not a random sample of the forest population: surveyors avoid inaccessible terrain and degraded stands, biasing observed AGB toward intermediate values. Standard multi-task and semi-supervised regression methods do not address this combination of structural label-partitioning, severe imbalance, and missingness-not-at-random (MNAR). We argue that large-scale forest attribute estimation is, at the level we attack here, a machine-learning problem with three coupled components. (i) The target variables must be predicted jointly, but their training labels come from disjoint sources covering disjoint subsets of the variable space, a condition stronger than standard partial-label learning. (ii) The plot labels are MNAR [Rubin, 1976]: a model fitted to the observed distribution will systematically underestimate AGB in dense forests. (iii) The targets are linked by an allometric relationship AGB ≈ g(H, C, SD, WD; ϕ) whose parameters ϕ vary by biome and cannot be fitted directly from any sample, since no sample contains all variables together. These characteristics should be treated jointly: balancing supervision sources in batches without correcting for MNAR amplifies the bias; correcting for MNAR without exploiting cross-variable structure leaves the unlabelled majority of pixels without supervision. We present StruMPL (Structured Multi-task Physics-constrained Learning), a single-objective framework that addresses the three problems jointly. A shared encoder produces features used by per-variable regression and imputation heads and by a propensity head that estimates per-pixel, per-variable observation probabilities. The regression objective is the Augmented Inverse Probability Weighting (AIPW) estimator [Robins et al., 1994, Scharfstein et al., 1999], with a stop-gradient design that lets the encoder be trained jointly with the propensity and imputation models without collapse. A differentiable physics module, with learnable allometric coefficients, evaluates the ecological law on the model’s own predictions at every pixel, providing semi-supervised gradient at the >99% of pixels that have no label. We evaluate the framework at two ecologically distinct sites: Mediterranean forests in Spain (SNFI [Ministerio para la Transición Ecológica y el Reto Demográfico, 2024]), and dry tropical/savanna forests in Africa (SEOSAW [SEOSAW-partnership, 2021]), demonstrating the framework’s portability across biomes and across differing variable availability. Contributions • We formalize multi-task dense regression under heterogeneous disjoint partial supervision with MNAR labels and inter-task physical constraints, a problem class appearing whenever structured domain knowledge links target variables observable only through separate, biased instruments. We test the formulation on AGB estimation across two biomes with different variable availability, using identical architecture and training recipe. • We introduce a joint propensity model for spatial MNAR correction in dense regression and adopt the AIPW pseudo-outcome estimator as the supervised loss, with separate imputation heads as the outcome baseline via stop-gradient. We show analytically and empirically that two stop-gradients (on the propensity and imputation baseline) are individually necessary for stable joint training: the propensity detach prevents propensity collapse, while the imputation 2
detach preserves the independence between outcome and imputation estimates that gives the AIPW pseudo-outcome its variance-reduction property over standard IPW. • We design a differentiable, parametrised constraint module that evaluates a known intertask biophysical relationship on the model’s own predictions at every pixel, supplying a cross-source semi-supervised signal at locations with no label. Its parameters are learned jointly with the network, and its output is initialised to a domain-plausible range to avoid the gradient pathologies that otherwise prevent joint training. We represent it as a learnable allometric form for forest biomass. • We validate the framework on two structurally distinct biomes (Mediterranean and dry tropical/savanna forests, with different target-variable availability) using identical architecture and training recipe. Stratified analysis shows the AIPW correction reduces bias specifically in the under-represented regime where MNAR effects concentrate; two distinct sites show that the unified mask design adapts cleanly when a target variable is missing.
2
Related Work
Forest attribute estimation from EO. AGB and canopy variables have been estimated from satellite imagery using radar empirical models [Ulaby et al., 1990] and, more recently, deep networks [Lang et al., 2023, Asiyabi et al., 2026, Weber et al., 2025]. Existing pipelines treat each variable as a separate target with independent training, or use disjoint per-source heads even when multiple variables are predicted [Asiyabi et al., 2026, Guo et al., 2023]. But none enforce the ecological constraint that links the variables, and none address the MNAR structure of plot labels. Multi-task learning under partial and heterogeneous supervision. Multi-task learning exploits inter-task statistical dependence to improve generalisation [Caruana, 1997], with task selection and gradient interference being central concerns [Standley et al., 2020, Yu et al., 2020]. Partiallabel MTL allows any subset of labels per sample [Nishi et al., 2024], and multi-source learning combines different supervision modalities [Ouyang et al., 2014, Hur et al., 2023, Fifty et al., 2021]. Both approaches assume a common label inventory across sources or random missingness; neither addresses the case where source identity determines which variables are labelled and the source distributions of the variables are linked by a known physical law. MNAR and propensity-based correction. The MCAR/MAR/MNAR are introduced in Rubin [1976] and subsequently developed by Little and Rubin [2019]. IPW [Horvitz and Thompson, 1952, Rosenbaum and Rubin, 1983] and AIPW [Robins et al., 1994, Scharfstein et al., 1999] are standard tools for unbiased estimation under MAR; the latter has a doubly-robust property under classical conditions. These estimators have been used in causal inference [Chernozhukov et al., 2018], recommendation [Schnabel et al., 2016, Saito et al., 2020], almost always at the sample level on tabular covariates. We apply AIPW at pixel resolution, with a per-pixel per-variable propensity estimated jointly with regression from a shared encoder, and characterise the gradient-separation needed for joint training. Physics-informed and semi-supervised learning. Physics-Informed Neural Networks [Raissi et al., 2019] encode physical laws as residual losses on a single target, typically with fixed parameters, and have been applied in fluid dynamics [Cai et al., 2021], climate [Kashinath et al., 2021], and ecology [Viet Cuong et al., 2024, Miranda et al., 2024]. Consistency regularisation [Tarvainen and Valpola, 2017, Sohn et al., 2020] provides another semi-supervised signal in dense prediction [Ouali et al., 2020, Ji et al., 2019]. Our physics module differs in two respects: its parameters are learned jointly with the network, so the model is not committed to published exponents that may not transfer; and it links disjoint label sources, providing a per-pixel pseudo-target derived from the model’s own structural predictions rather than from a teacher network or a confident-prediction threshold.
3
Problem Formulation
Setup. Let x ∈ RCin ×H×W be a multi-channel EO patch, y ∈ RK×H×W a joint target tensor over K forest attributes, and R ∈ {0, 1}K×H×W a binary observation mask with Rk,i,j = 1 if variable k 3
is labelled at pixel (i, j). The goal is a mapping fθ : x 7→ ŷ that minimises expected error over all variables and pixels in the target population. Heterogeneous disjoint and sparse label structure. Training data is sparse (only a few labelled pixels in each patch) and comes from two sources. Spaceborne lidar (GEDI) provides KG ⊂ {1, . . . , K} at sparse ∼25 m footprints, with NG samples; field plots provide KP ⊂ {1, . . . , K} at plot locations, with NP samples. By construction KG ∩ KP = ∅: GEDI carries {H, C} and plots carry {SD, AGB, WD}; the source determines which variables are labelled. The full training set is D = DG ∪ DP , with severe imbalance NG ≫ NP (∼250:1 for Africa, ∼110:1 for Spain). This is stronger than standard partial-label learning, where the missingness is not deterministic by the source. MNAR labels. Within DP , the labelled locations are not a random sample. Surveyors visit accessible forests at moderate elevation and avoid both inaccessible old-growth and degraded stands; the high-AGB tail is therefore underrepresented. Following Rubin [1976], labels are Missing Not At Random (MNAR) when P (Rk,i,j = 1 | x, y) depends on yk,i,j even given x. A model trained by the standard naive masked supervised loss optimises performance on the observed distribution rather than the population, and in our setting will systematically underestimate dense forests and overestimate degraded ones. The standard correction is to model the propensity πk,i,j = P (Rk,i,j = 1 | x) and re-weight by 1/πk,i,j [Rosenbaum and Rubin, 1983]. This recovers an unbiased estimator under conditional ignorability (i.e., that label availability depends on x but not on yk,i,j given x). In our setting, however, ignorability does not hold exactly as the physical inaccessibility may correlate with biomass beyond what the EO data show at the patch scale. Our assumption is that the encoder’s learned representation is a sufficiently rich covariate summary for ignorability to hold approximately. To reduce sensitivity to propensity misspecification we use the AIPW estimator [Robins et al., 1994], which has a doubly-robust property under our working assumption of approximate ignorability: it is consistent if either the propensity or the imputation model is correctly specified. Physical constraint.
Forest ecology gives an allometric relationship AGB ≈ g H, C, SD, WD; ϕ , (1) with biome-specific parameters ϕ. Crucially, equation (1) cannot be fitted from the training data alone as KG ∩ KP = ∅ (i.e., H, C are never co-observed with AGB, SD, WD on any sample). We therefore introduce g as an external structural prior. Once fθ produces ŷ, the allometric residual ŷAGB − g(ŷH , ŷC , ŷSD , ŷWD ; ϕ) is defined at every pixel, labelled or not. This makes the physical constraint a natural cross-source semi-supervised signal at the >99% of pixels with no label. Multi-site variable availability. Survey protocols differ in which attributes they record: SEOSAW (Africa) reports KP = {SD, WD, AGB} while SNFI (Spain) lacks wood density (KP = {SD, AGB}). Variable absence is encoded as an all-zero mask column, so the architecture is unchanged across sites; only the input set of the physics function g changes. The sites also differ in test-set size: SNFI provides ∼6,000 plots, supporting stratified bias analysis; SEOSAW provides 162 standardised 1 ha plots recommended by the data providers as the network’s higher-quality subset (Appendix C and K). Stratified analysis is therefore reported only on SNFI; aggregate metrics are reported on both.
4
Method
StruMPL addresses the three problems above within a single forward pass and a unified objective (Figure 1). A shared encoder Encθ : x 7→ z ∈ RD×H×W produces features used by K independent regression heads (each producing ŷk ∈ RH×W in z-score space), K independent imputation heads (each producing m̂k ∈ RH×W in z-score space) and a propensity head π̂ = σ(BiasHead(z)) ∈ (0, 1)K×H×W . We use a ResUNet with channel and spatial attention as the encoder; full architectural details are in Appendix A. The total objective is Ltotal = Lsup + λphys Lphys + λcons Lcons + λbias Lbias + λimp Limp , (2) with each term targeting one of the failure modes identified in section 3. We describe the two terms central to the contribution (Lsup and Lphys ) below; Lcons , Lbias and Limp are standard and we summarise them in section 4.3 (full forms in Appendix B). 4
Figure 1: Overview of StruMPL. A shared encoder feeds regression, imputation and propensity heads. Lsup uses the AIPW pseudo-outcome with stop-gradients on π̂ and m̂ (section 4.1); Limp trains the imputation heads on the labelled pixels; Lphys enforces allometric consistency at every pixel via the learnable module gϕ ; Lcons regularises under input augmentation; Lbias trains π̂ against the mask. 4.1
Supervised loss with AIPW MNAR correction
The naive masked MSE optimises performance on the observed distribution. We replace it with the AIPW pseudo-outcome estimator. Let µk,i,j = sg(m̂k,i,j ) denote the imputation estimate with gradients stopped, and π̂k,i,j the propensity (also gradients stopped), clamped to [πmin , 1] with πmin = 0.1 to bound the effective AIPW weight. We define the pseudo-outcome ỹk,i,j = µk,i,j +
Rk,i,j yk,i,j − µk,i,j . π̂k,i,j
(3)
The supervised loss is a pixel-wise MSE between the regression output and the pseudo-outcome: X 2 1 Lsup = ŷk,i,j − ỹk,i,j . (4) K · HW k,i,j
Stop-gradients are not optional. Two stop-gradients in Eq. (3) keep joint training stable. Without detaching π̂, the encoder could minimise Lsup by inflating π̂ → 1, collapsing the IPW weight and erasing the MNAR correction. Without detaching µ, Lsup would backpropagate into the imputation heads and drive m̂ → ŷ, collapsing Lsup to a 1/π̂ 2 -weighted MSE on labelled pixels (worse than standard IPW). The propensity and imputation heads are therefore trained only through Lbias and Limp (section 4.3), and enter Lsup as fixed coefficients. With both stop-gradients applied, the gradient w.r.t. ŷ is 2(ŷ − ỹ), with no gradient flowing into m̂ or π̂. At unlabelled pixels (R = 0), this pulls ŷ toward µ; at labelled pixels (R = 1), the pseudo-outcome combines µ with an IPW-corrected residual, so the model is pushed toward the AIPW-weighted stationary points. When µ is reasonably well-fit, the residual (y − µ) is small, reducing the O(1/π̂) loss spikes characteristic of standard IPW; this is the variance-reduction property of AIPW over IPW (full derivation in Appendix E). 5
We note one caveat: classical double robustness requires cross-fitting for asymptotic independence of the estimators [Chernozhukov et al., 2018]. Due to our shared-encoder architecture, we use AIPW for empirical robustness to either-side misspecification rather than as a formal consistency guarantee. 4.2
Physics consistency loss
The physics module gϕ is a learnable allometric equation evaluated on the model’s own structural predictions. We use the ecology-motivated form [ phys = α + scale · sp(ŷH )b · sp(ŷC )c sp(ŷSD )·d · sp(ŷWD )e , AGB
(5)
where ϕ = {α, scale, b, c, d, e} are constrained positive via softplus (sp(·) = log(1 + exp(·))) or exp. The interaction term [sp(ŷH )b · sp(ŷC )c ]sp(ŷSD )·d allows stem density to modulate the structural-tobiomass scaling, motivated by competition effects in dense stands. When WD is unavailable (Spain), the sp(ŷWD )e factor is dropped and |ϕ| decreases by one. The form is an ecologically-motivated parametrisation rather than a derivation; we compare it against a plain multiplicative power law and a learned MLP in section 5. The regression heads operate in z-score space; the physics module operates in physical units. We therefore denormalise the structural predictions before equation (5), clamp the output to the ecologically plausible range [0, 2000] Mg/ha, and renormalise to z-score space using the training-set statistics. The physics consistency loss is then a pixelwise MSE applied to all pixels, labelled or not: Lphys =
norm 1 X [ phys,i,j 2 . ŷAGB,i,j − AGB HW i,j
(6)
Gradients flow from Lphys back into all structural heads, into the AGB head, and into ϕ. This is the mechanism by which the two label sources exchange information: GEDI-supervised ŷH , ŷC are constrained, under the learned allometry, to be consistent with plot-supervised ŷAGB . 4.3
Other components
Propensity training. The BiasHead is trained via binary cross-entropy between π̂ and the observation mask R, applied at every pixel. Gradients flow back through the shared encoder, so the encoder is trained to represent features predictive of survey accessibility, i.e., exactly the features the propensity model needs. The full form is in Appendix B. Augmentation consistency. We add a standard FixMatch-style consistency term [Tarvainen and Valpola, 2017, Sohn et al., 2020]: different predictions on x and on a perturbed view x′ (additive Gaussian noise and per-element channel dropout) are penalised. We do not use EMA teacher and confidence threshold, as the small patch size and the stronger physics signal make these unnecessary. Imputation training. The imputation heads m̂k are trained via masked MSE on labelled pixels only, providing the outcome baseline used in the AIPW pseudo-outcome (section 4.1). They are deliberately trained on the (biased) observed distribution: this is what AIPW expects of the outcome model, and the propensity correction in Lsup debiases the final regression heads. Implementation details. Complete implementation details including architecture components, hyperparameters, source-balanced batching, parameters initialisation, and parameter count are given in Appendix A, C, and D.
5
Results and Discussion
We evaluate StruMPL on two datasets covering ecologically distinct biomes: Spain (SNFI, ∼6,000 test plots) and Africa (SEOSAW, 162 test plots in the standardised 1 ha format). For headline numbers we report mean and standard deviation over 5 random seeds (42, 123, 456, 789, 1011); ablation rows use a single seed (42) to keep compute tractable. 6
Table 1: Main results on the primary AGB target (Mg/ha). Lower is better for both RMSE and |bias|. Full StruMPL reports mean ± std over 5 seeds; ablations are single-seed. Spain (SNFI) Configuration
Africa (SEOSAW)
AGB RMSE
AGB bias
AGB RMSE
AGB bias
External baseline ([Santoro and Cartus, 2024]) External baseline ([Asiyabi et al., 2026]) Single-output AGB, naive MSE Single-output AGB + Lbias + Lcons + Limp
65.3 44.8 43.7 43.1
18.8 1.2 0.9 2.1
38.5 21.8 22.5 22.6
27.9 1.5 5.3 3.0
Joint multi-task variants MTL, sup only (no Lphys , Lcons , Lbias , Limp ) MTL + Lphys only MTL + Lbias + Limp MTL + Lcons only Full StruMPL (AIPW + allometric)
43.9 43.7 43.0 43.3 41.9 ± 0.5
1.3 0.9 1.5 0.8 0.4 ± 0.1
24.9 24.1 23.6 23.7 20.7 ± 0.6
1.4 5.4 2.2 1.3 0.6 ± 0.7
Figure 2: AGB RMSE (bars, left axis) and AGB bias (red markers, right axis) for all model configurations from Table 1 on the Spain (left) and Africa (right) test sets. Lower is better for both metrics. Error bars on Full StruMPL denote standard deviation across 5 random seeds. 5.1
Headline performance and the cost of each component
Table 1 and Figure 2 report AGB RMSE and bias for full StruMPL, MTL variants with loss-component ablations, two single-output baselines, and two published external baselines ([Asiyabi et al., 2026] closest to StruMPL, and [Santoro and Cartus, 2024] standard global AGB map from ESA CCI; see Appendix F). The complete ablation grid and per-variable RMSE for the structural variables are in Appendix I; here we focus on AGB as the main target variable. Discussion. Full StruMPL achieves the best AGB RMSE and bias on both sites, but the margins tell different stories. On Africa, where the source imbalance is larger, StruMPL shows a clear gain over every ablation: 20.7 vs. 22.5 Mg/ha against the strongest single-component variant, and 24.9 Mg/ha against the MTL baseline with no propensity correction or physics consistency. The MTL-only configuration is in fact the worst row on Africa, showing that MTL without the StruMPL components is actively harmful, and that the observed gains are not a generic side-effect of multi-tasking but specifically attributable to Lphys and Lbias . The Africa results also support the portability claim: the unified mask design adapts cleanly to wood density availability, the BalancedBatchSampler handles the larger source imbalance (NG /NP ≈ 250 vs. Spain’s ≈ 110), and the curriculum warmup is effective without recipe changes (further discussed in Appendix K). On Spain, the gain over the strongest ablation is smaller but statistically robust: full StruMPL reaches 41.9 ± 0.5 Mg/ha against 43.0–43.9 for the MTL variants and 43.1 for the strongest single-output baseline. The improvement over the strongest ablation (MTL + Lbias + Limp ) is significant under a paired bootstrap (p < 0.05, n = 10,000; Appendix H). The clearer StruMPL signal on Spain is in bias: full StruMPL achieves |bias| = 0.4 ± 0.1 Mg/ha against 0.8–2.1 for the ablations, which is the impact of the AIPW correction (section 5.2). StruMPL also outperforms the closest published method [Asiyabi et al., 2026] on both sites in both metrics, and is substantially more accurate than 7
the off-the-shelf ESA CCI global biomass map [Santoro and Cartus, 2024], which we include as a reference point. 5.2
Where AIPW pays off: stratified bias on Spain
The AIPW correction is designed to reduce systematic high-AGB underestimation where MNAR effects concentrate. Aggregate metrics hide this because the high-AGB tail is a small fraction of test samples. We bin Spain test pixels (due to the higher number of the test plots) into 5 AGB quantiles by plot label and report RMSE and bias per quantile (Figure 3).
Figure 3: Bias and RMSE on AGB across five quantiles for naive masked MSE, IPW, and AIPW supervision (all other components fixed). AIPW reduces the systematic high-AGB underestimation that naive MSE exhibits; per-quantile RMSE redistributed toward the under-represented high tail. Discussion. The naive masked MSE shows the predicted MNAR pattern: overestimation in low AGB regions (+25 Mg/ha in Q1) and large underestimation in the high-AGB tail (−50 Mg/ha in Q5), crossing zero around Q4; whole-dataset bias (+1.7 Mg/ha) hides this entirely as the two ends cancel. IPW reduces Q5 underestimation modestly (−40.6 Mg/ha) while AIPW reduces it substantially (−22.8 Mg/ha, a 54% reduction); lower quantiles improve proportionally (Q1: +25 → +10.1). The RMSE pattern reflects the same redistribution: AIPW reduces Q5 RMSE from 65.4 to 49.2 Mg/ha while marginally increasing RMSE in Q1–Q4 (by 3–4 Mg/ha), as expected from reweighting toward under-represented samples. The net effect on whole-dataset RMSE is positive (43.2 to 41.9 Mg/ha), and the framework prioritises the regime where MNAR error is concentrated. 5.3
Are the stop-gradients necessary?
In section 4.1 and Appendix E, we mentioned that two stop-gradients in the AIPW pseudo-outcome are individually necessary for stable joint training. We test this empirically (Table 2). Table 2: Stop-gradient ablation on the Spain test set, full StruMPL configuration. Configuration Both stop-gradients (default) π̂ not detached in Lsup µ = ŷ not detached
AGB RMSE
Mean π̂ at convergence
Outcome
41.9 44.2 44.6
0.62 0.99 0.60
stable propensity collapse imputation heads collapse
Discussion. Without detaching π̂, the encoder reduces Lsup by artificially inflating predicted propensities; we observe mean π̂ drift to ≈ 1.0, which collapses the AIPW correction into a naive, biased MSE. Without detaching µ, the imputation heads lose their functional independence and begin to shadow the main regression heads. This results in the 1/π̂ 2 weighting failure described in Appendix E, which amplifies noise and degrades performance. These empirical failures confirm that the stop-gradients are structural requirements of the AIPW framework. 5.4
Ecologically motivated physics module
In Table 3, we compare the chosen allometric parametrisation against power_law (a plain multiplicative power law, the conservative baseline) and mlp (a 3-layer MLP, an unconstrained baseline). 8
Table 3: Physics module form; full StruMPL otherwise. Physics module
Spain AGB RMSE
Africa AGB RMSE
Notes
44.3 41.9 44.4
25.5 20.7 24.2
baseline ecology motivated not interpretable
power_law allometric (default) mlp
Discussion. The allometric form outperforms both alternatives on both sites, reducing AGB RMSE by 2.4 Mg/ha (5.4%) over power_law on Spain and 4.8 Mg/ha (18.8%) on Africa. A 3layer MLP has more parameters and more flexibility than either parametrised form, yet performs comparably to the plain power_law on Spain (44.4 vs. 44.3) and remains worse than allometric on Africa. The improvement therefore comes from inductive bias rather than expressiveness. The relative gain is larger on Africa than on Spain. This observation is consistent with the ecological intuition that competition-driven scaling effects matter more in the heterogeneous-density savanna stands of the SEOSAW domain than in the more homogeneous Mediterranean stands of SNFI; we report this as a qualitative observation rather than a quantitative claim, given the smaller Africa test set. We use allometric as the default and use its learned coefficients in the cross-site comparison of section 5.5. 5.5
Learned allometric coefficients
The allometric parameters ϕ = {α, scale, b, c, d, e} are jointly trained via Lphys . Their fitted values are not maximum-likelihood ecological estimates as they are optimised under the full StruMPL objective (including AIPW reweighting and consistency regularisation), but they characterise the effective satellite-scale scaling that the model converges to. Table 4 reports them across sites. Table 4: Learned allometric exponents. Spain has no WD term (e absent). Cross-site comparison is qualitative; see text. Site
α
scale
b (H)
c (C)
d (SD)
e (WD)
Spain Africa
0.82 0.54
31.3 21.0
0.04 0.05
0.21 0.03
0.02 0.03
— 0.02
Discussion. Several points stand out in Table 4. First, the structural exponents (b, c, d, e) are uniformly small, far from the near-unit values of textbook allometric equations [Chave et al., 2014]. We do not interpret this as a recovery of the textbook allometry: under the StruMPL objective, the optimiser finds a regime in which AGB is largely explained by the additive baseline α and the multiplicative scale, with the structural exponents acting as small corrections. This matches the framing of section 4.2: ϕ are effective satellite-scale scaling coefficients, not estimates of the underlying ecological law. Second, the cross-site differences are concentrated in two parameters. The Spain scale (31.3) exceeds Africa’s (21.0), consistent with the higher mean AGB density of Mediterranean stands; and the canopy-cover exponent c is markedly higher on Spain (0.21 vs. 0.03), aligning with the intuition that denser-canopy Mediterranean forests show a stronger biomasscover relationship than savanna stands, where local cover can be high (an isolated tree) without corresponding landscape biomass. The remaining exponents are roughly comparable across sites. We treat these readings as qualitative ecological observations rather than quantitative estimates. 5.6
Limitations
StruMPL relies on two assumptions. First, conditional ignorability (i.e., the labels are MAR given the encoder’s learned representation, i.e., P (R = 1 | z, y) = P (R = 1 | z)), however, in our setting the labels are MNAR with respect to the raw inputs x (plot accessibility correlates with biomass), and we rely on encoder to absorb enough of this dependence into z that conditional independence holds approximately. AIPW does not protect against violations of this assumption itself; it protects only against misspecification of the propensity or the imputation model when ignorability holds. Second, the formal double-robustness guarantee for ML estimators requires cross-fitting to control overfitting bias [Chernozhukov et al., 2018]. We share the encoder between the heads and do not cross-fit, so 9
we use AIPW as an empirical bias-reduction tool rather than a formal consistency guarantee. Both points constrain the theoretical claims we can make but do not invalidate the empirical bias reduction reported in section 5.2. A separate caveat is the choice of allometric form: ours is one of several plausible parametrisations, and the learned exponents should be read as effective scaling coefficients under the StruMPL objective, not ecological estimates. We compare against alternative forms in section 5.4.
6
Conclusion
We presented StruMPL, a framework for multi-task dense regression under disjoint heterogeneous partial supervision with MNAR labels and inter-task physical constraints. The framework combines a joint propensity model for spatial MNAR correction, an AIPW pseudo-outcome loss with two stop-gradients whose necessity we establish both analytically and empirically, and a learnable physics module that provides cross-source semi-supervised gradient at every pixel. Applied to forest attribute estimation on two ecologically distinct biomes, StruMPL outperforms single-output baselines, ablation variants, and the closest published method, with the AIPW correction specifically reducing bias in the high-AGB regime where MNAR effects are concentrated. The four-condition characterisation of when the formulation applies (Appendix M) suggests that domains with structurally similar properties (hydrology, battery state estimation, multi-modal medical imaging) are candidate applications. We leave cross-domain test and addressing the limitations (section 5.6) to future work.
10
References Reza M Asiyabi, SEOSAW Partnership, Steven Hancock, and Casey Ryan. Process-guided concept bottleneck model. arXiv preprint arXiv:2601.10562, 2026. Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Physicsinformed neural networks (PINNs) for fluid mechanics: A review. Acta Mechanica Sinica, 37(12): 1727–1738, 2021. Rich Caruana. Multitask learning. Machine learning, 28(1):41–75, 1997. Jérôme Chave, Maxime Réjou-Méchain, Alberto Búrquez, Emmanuel Chidumayo, Matthew S Colgan, Welington BC Delitti, Alvaro Duque, Tron Eid, Philip M Fearnside, Rosa C Goodman, et al. Improved allometric models to estimate the aboveground biomass of tropical trees. Global change biology, 20(10):3177–3190, 2014. Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, Whitney Newey, and James Robins. Double/debiased machine learning for treatment and structural parameters, 2018. Ralph Dubayah, James Bryan Blair, Scott Goetz, Lola Fatoyinbo, Matthew Hansen, Sean Healey, Michelle Hofton, George Hurtt, James Kellner, Scott Luthcke, et al. The global ecosystem dynamics investigation: High-resolution laser ranging of the earth’s forests and topography. Science of remote sensing, 1:100002, 2020. Ted R Feldpausch, Jon Lloyd, Simon L Lewis, Roel JW Brienen, Manuel Gloor, Abel Monteagudo Mendoza, Gabriela Lopez-Gonzalez, Lindsay Banin, Kamariah Abu Salim, Kofi AffumBaffoe, et al. Tree height integrated into pantropical forest biomass estimates. Biogeosciences, 9 (8):3381–3403, 2012. Chris Fifty, Ehsan Amid, Zhe Zhao, Tianhe Yu, Rohan Anil, and Chelsea Finn. Efficiently identifying task groupings for multi-task learning. Advances in Neural Information Processing Systems, 34: 27503–27516, 2021. Qiyu Guo, Shouhang Du, Jinbao Jiang, Wei Guo, Hengqian Zhao, Xuzhe Yan, Yinpeng Zhao, and Wanshan Xiao. Combining GEDI and sentinel data to estimate forest canopy mean height and aboveground biomass. Ecological Informatics, 78:102348, 2023. Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016. Daniel G Horvitz and Donovan J Thompson. A generalization of sampling without replacement from a finite universe. Journal of the American statistical Association, 47(260):663–685, 1952. Jie Hu, Li Shen, and Gang Sun. Squeeze-and-excitation networks. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 7132–7141, 2018. Kyunghoon Hur, Jungwoo Oh, Junu Kim, Jiyoun Kim, Min Jae Lee, Eunbyeol Cho, Seong-Eun Moon, Young-Hak Kim, Louis Atallah, and Edward Choi. Genhpf: General healthcare predictive framework for multi-task multi-source learning. IEEE Journal of Biomedical and Health Informatics, 28(1):502–513, 2023. IPCC Task Force on National Greenhouse Gas Inventories (TFI). 2019 refinement to the 2006 IPCC guidelines for national greenhouse gas inventories. https://www.ipcc-nggip.iges.or.jp/ public/2019rf/index.html, 2019. Accessed: 2026-04-15. Japan Aerospace Exploration Agency (JAXA). ALOS PALSAR-2 ScanSAR Level 2.2 Dataset. https://developers.google.com/earth-engine/datasets/catalog/JAXA_ ALOS_PALSAR-2_Level2_2_ScanSAR, 2021. Provided by EORC, accessed through Google Earth Engine. Rongrong Ji, Ke Li, Yan Wang, Xiaoshuai Sun, Feng Guo, Xiaowei Guo, Yongjian Wu, Feiyue Huang, and Jiebo Luo. Semi-supervised adversarial monocular depth estimation. IEEE transactions on pattern analysis and machine intelligence, 42(10):2410–2422, 2019. 11
Karthik Kashinath, Mustafa Mustafa, Adrian Albert, Jean-Luc Wu, C Jiang, Soheil Esmaeilzadeh, Kamyar Azizzadenesheli, R Wang, Ashesh Chattopadhyay, Aakanksha Singh, et al. Physicsinformed machine learning: case studies for weather and climate modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 379(2194), 2021. Nico Lang, Walter Jetz, Konrad Schindler, and Jan Dirk Wegner. A high-resolution canopy height model of the earth. Nature Ecology & Evolution, 7(11):1778–1789, 2023. Roderick JA Little and Donald B Rubin. Statistical analysis with missing data. John Wiley & Sons, 2019. Ilya Loshchilov and Frank Hutter. arXiv:1711.05101, 2017.
Decoupled weight decay regularization.
arXiv preprint
Jeffrey Masek, Junchang Ju, Jean-Claude Roger, Sergii Skakun, Eric Vermote, Martin Claverie, Jennifer Dungan, Zhangshi Yin, Brian Freitag, and Chris Justice. HLS operational land imager surface reflectance and TOA brightness daily global 30m v2. 0. NASA EOSDIS Land Processes Distributed Active Archive Center (DAAC) data set, pages HLSL30–002, 2021. Ministerio para la Transición Ecológica y el Reto Demográfico. Cuarto Inventario Forestal Nacional (IFN4). Ministerio para la Transición Ecológica y el Reto Demográfico, 2024. ISBN 978-84-18778-51-3. URL https://libreria.miteco.gob.es/libro/ cuarto-inventario-forestal-nacional_4635/. Miro Miranda, Marcela Charfuelan, and Andreas Dengel. Exploring physics-informed neural networks for crop yield loss forecasting. arXiv preprint arXiv:2501.00502, 2024. Kento Nishi, Junsik Kim, Wanhua Li, and Hanspeter Pfister. Joint-task regularization for partially labeled multi-task learning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 16152–16162, 2024. Ozan Oktay, Jo Schlemper, Loic Le Folgoc, Matthew Lee, Mattias Heinrich, Kazunari Misawa, Kensaku Mori, Steven McDonagh, Nils Y Hammerla, Bernhard Kainz, et al. Attention u-net: Learning where to look for the pancreas. arXiv preprint arXiv:1804.03999, 2018. Yassine Ouali, Céline Hudelot, and Myriam Tami. Semi-supervised semantic segmentation with cross-consistency training. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 12674–12684, 2020. Wanli Ouyang, Xiao Chu, and Xiaogang Wang. Multi-source deep learning for human pose estimation. In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2329– 2336, 2014. Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. PyTorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019. Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019. Daniela Requena Suarez, Danaë MA Rozendaal, Veronique De Sy, Oliver L Phillips, Esteban AlvarezDávila, Kristina Anderson-Teixeira, Alejandro Araujo-Murakami, Luzmila Arroyo, Timothy R Baker, Frans Bongers, et al. Estimating aboveground net biomass change for tropical and subtropical forests: Refinement of IPCC default rates using forest plot data. Global Change Biology, 25 (11):3609–3624, 2019. James M Robins, Andrea Rotnitzky, and Lue Ping Zhao. Estimation of regression coefficients when some regressors are not always observed. Journal of the American statistical Association, 89(427): 846–866, 1994. 12
Paul R Rosenbaum and Donald B Rubin. The central role of the propensity score in observational studies for causal effects. Biometrika, 70(1):41–55, 1983. Donald B Rubin. Inference and missing data. Biometrika, 63(3):581–592, 1976. Yuta Saito, Suguru Yaginuma, Yuta Nishino, Hayato Sakata, and Kazuhide Nakata. Unbiased recommender learning from missing-not-at-random implicit feedback. In Proceedings of the 13th international conference on web search and data mining, pages 501–509, 2020. Mattia Santoro and Oliver Cartus. ESA Biomass Climate Change Initiative (Biomass CCI): Global datasets of forest above-ground biomass for the years 2010–2021, v5.01, 2024. URL https: //doi.org/10.5285/bf535053562141c6bb7ad831f5998d77. Accessed: 2026-02-05. Daniel O Scharfstein, Andrea Rotnitzky, and James M Robins. Adjusting for nonignorable drop-out using semiparametric nonresponse models. Journal of the American Statistical Association, 94 (448):1096–1120, 1999. Tobias Schnabel, Adith Swaminathan, Ashudeep Singh, Navin Chandak, and Thorsten Joachims. Recommendations as treatments: Debiasing learning and evaluation. In international conference on machine learning, pages 1670–1679. PMLR, 2016. SEOSAW-partnership. A network to understand the changing socio-ecology of the southern african woodlands (SEOSAW): Challenges, benefits, and methods. Plants, People, Planet, 3(3):249–267, 2021. Kihyuk Sohn, David Berthelot, Nicholas Carlini, Zizhao Zhang, Han Zhang, Colin A Raffel, Ekin Dogus Cubuk, Alexey Kurakin, and Chun-Liang Li. Fixmatch: Simplifying semi-supervised learning with consistency and confidence. Advances in neural information processing systems, 33:596–608, 2020. Trevor Standley, Amir Zamir, Dawn Chen, Leonidas Guibas, Jitendra Malik, and Silvio Savarese. Which tasks should be learned together in multi-task learning? In International conference on machine learning, pages 9120–9132. PMLR, 2020. Antti Tarvainen and Harri Valpola. Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results. Advances in neural information processing systems, 30, 2017. Fawwaz T Ulaby, Kamal Sarabandi, K.Y.L.E. Mcdonald, Michael Whitt, and M Craig Dobson. Michigan microwave canopy scattering model. International Journal of Remote Sensing, 11(7): 1223–1253, 1990. Dinh Viet Cuong, Branislava Lalić, Mina Petrić, Nguyen Thanh Binh, and Mark Roantree. Adapting physics-informed neural networks to improve ODE optimization in mosquito population dynamics. PLOS One, 19(12):e0315762, 2024. Manuel Weber, Carly Beneke, and Clyde Wheeler. Unified deep learning model for global prediction of aboveground biomass, canopy height, and cover from high-resolution, multi-sensor satellite imagery. Remote Sensing, 17(9):1594, 2025. Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn. Gradient surgery for multi-task learning. Advances in neural information processing systems, 33: 5824–5836, 2020.
13
Appendix A
Architecture
Shared encoder. The encoder Encθ is a ResUNet [He et al., 2016] with residual blocks, Squeezeand-Excitation channel attention [Hu et al., 2018], and spatial attention gates on the decoder skip connections [Oktay et al., 2018]. Given x ∈ RCin ×H×W with H = W = 16, it produces z = Encθ (x) ∈ RD×H×W with D = 128. The encoder follows the standard U-Net contractive–expansive structure with 3 encoder stages, channel widths doubling at each downsampling level, and matching upsampling in the decoder. The shared encoder is motivated by two observations specific to this problem. First, NG ≫ NP : borrowing statistical strength from the larger GEDI corpus is only possible under a shared representation. Second, the physics consistency loss links variables supervised by different sources, and the residual is meaningful only if the structural and biomass predictions are computed from a compatible feature space. The input contains Cin = 15 channels drawn from a multi-sensor EO time series: ALOS PALSAR2 L-band SAR backscatter [Japan Aerospace Exploration Agency (JAXA), 2021], Harmonized Landsat Sentinel-2 (HLS) optical reflectance bands [Masek et al., 2021], and geographic coordinates. Channel attention lets the encoder learn input-channel importance per task; spatial attention on skip connections suppresses uninformative spatial regions, useful given that labels cover <1% of pixels per patch. Regression heads.
K independent two-layer convolutional heads map z to per-pixel predictions, k ∈ {1, . . . , K}.
ŷk = Headk (z),
(7)
Each head consists of a 3 × 3 convolution (D → D/2 channels) with RELU activation, followed by a 1 × 1 convolution (D/2 → 1). We keep heads independent because each variable receives gradient from a different mask and a different source; sharing parameters would induce gradient interference between tasks with incompatible observability. Inter-variable correlation is enforced through Lphys , which provides cross-variable signal without requiring shared head parameters. All predictions are produced in z-score normalised space (zero mean, unit variance per variable, computed on the training set). This makes Lsup scale-invariant across variables with different physical ranges, makes Lphys numerically comparable to Lsup , and keeps gradient magnitudes in a predictable range. Imputation heads. Similar to the regression heads, K independent two-layer convolutional heads with the same architecture (3 × 3 conv D → D/2 with RELU, followed by 1 × 1 conv D/2 → 1) produce the per-pixel imputation outputs, m̂k = ImputationHeadk (z),
k ∈ {1, . . . , K}.
(8)
The imputation heads share the encoder with the regression and bias heads, but are kept architecturally separate from the regression heads, with no parameter sharing between ŷk and m̂k . This separation preserves the functional independence between the outcome and imputation estimators that gives AIPW its variance-reduction property over standard IPW (section 4.1). Like the regression heads, imputation outputs are produced in z-score normalised space using the same training-set statistics, so µk = sg(m̂k ) enters Lsup on the same scale as ŷk and yk . The imputation heads are trained exclusively through Limp (Appendix B); gradients from Lsup are stopped at m̂k , preventing the imputation model from collapsing onto the regression head (Appendix E). Bias head.
The BiasHead is a small convolutional head with sigmoid output, π̂ = σ(BiasHead(z)) ∈ (0, 1)K×H×W .
(9)
Its architecture mirrors the regression and imputation heads (3 × 3 conv D → D/2 with RELU, 1 × 1 conv D/2 → K) so that the propensity estimation capacity matches the regression capacity. It shares the encoder with the regression heads, so the features that drive prediction also drive propensity estimation. 14
Parameter count. The full StruMPL model has ∼16.9 M parameters, of which the encoder accounts for 16,204,198, the regression heads for 369,605 (73,921 each head), the imputation heads also for 369,605 (73,921 each head), the BiasHead for 74,181, and the physics module for 6 (a constant; the physics parameters are scalar). Removing the imputation heads, BiasHead and physics module gives a standard multi-task ResUNet.
B
Full forms of Lbias , Lcons and Limp
Propensity training loss. Lbias = −
X 1 Rk,i,j log π̂k,i,j + (1 − Rk,i,j ) log(1 − π̂k,i,j ) . K · HW
(10)
k,i,j
This loss is evaluated at every pixel: labelled pixels are positive examples, unlabelled pixels are negative. Gradients flow back through the BiasHead into the encoder, training it to represent features correlated with label availability, exactly the features that the propensity model needs for downstream MNAR correction. Augmentation consistency loss. 1 X 2 ŷi,j (x) − ŷi,j (x′ ) , Lcons = HW i,j
x′ = x ⊙ B + ε,
(11)
with ε ∼ N (0, σ 2 I) (σ = 0.05) and B ∼ Bernoulli(1 − pdrop )Cin ×H×W (pdrop = 0.05). The perturbations are semantically neutral for EO patches at 16×16 resolution. We do not use Exponential Moving Average (EMA) teacher and confidence thresholding (cf. FixMatch [Sohn et al., 2020], Mean Teacher [Tarvainen and Valpola, 2017]); these refinements are standard in semi-supervised image classification but unnecessary here, given the small patch size (16×16) and the stronger cross-variable signal from Lphys . The augmented forward pass skips the imputation heads, BiasHead and physics module since they do not contribute to Lcons . Imputation training loss. X 2 1 Rk,i,j m̂k,i,j − yk,i,j . R k,i,j k,i,j
Limp = P
(12)
k,i,j
This is naive masked MSE between the imputation head outputs m̂k,i,j and the labels yk,i,j , applied only at labelled pixels. The imputation heads are trained exclusively through this loss: gradients from Lsup are stopped at m̂ (cf. equation (3)), so the outcome estimator and the regression head do not share gradient flow. This preserves the independence between the outcome model and the propensity correction that gives AIPW its variance-reduction property over standard IPW. The imputation heads are themselves biased toward the observed distribution (since they are trained on labelled pixels only); this is the standard AIPW design, and the bias is corrected by the propensity weighting in Lsup . The architecture of the imputation heads matches the regression heads (Appendix A), giving the imputation model the same representational capacity as the deployed regression model.
C
Source-balanced sampling
Without intervention, a uniform sampler over D = DG ∪ DP would draw samples in proportion to source size: ∼99.6% GEDI / 0.4% plot for Africa, ∼99.1% / 0.9% for Spain. The plot-supervised heads (AGB, SD, WD) would then receive gradient updates in fewer than 1% of batches, and the physics consistency loss, which requires GEDI-calibrated H, C and plot-calibrated AGB to coexist in a single batch, would rarely receive a balanced cross-source pairing. As a result, we use a source-balanced sampler where every batch is constructed as a 50/50 mix: G P P Batch = xG 1 , . . . , xB/2 ∪ x1 , . . . , xB/2 . | | {z } {z } GEDI
(13)
Plot
GEDI samples are shuffled without replacement per epoch; plot samples are oversampled with replacement. Epoch length is ⌊NG /(B/2)⌋, so each GEDI sample is seen approximately once per epoch and each plot sample is seen NG /NP times. 15
Effective passes per epoch. Because plot data are heavily oversampled, the two sources converge on different timescales:
Training samples Test samples Effective passes per epoch Physics warmup bottleneck?
GEDI (H, C)
Plot, Spain
Plot, Africa
∼2,000,000 ∼10,000 1 yes
∼20,000 ∼6,000 ∼110 no
∼8,000 162 ∼250 no
The physics warmup is therefore governed by epoch count (which tracks GEDI convergence rather than step count or plot pass count). Test set size. The test set sizes (Spain: ∼6,000 plots covering the region of interest; Africa: 162 standardised 1 ha plots, recommended by the SEOSAW data owners as the standard and higher-quality plots) further motivate the asymmetric experimental treatment in section 5: stratified analyses are reported on Spain, while Africa supports aggregate-level claims about portability. Why 50/50. Lower plot ratios under-represent the only source of AGB supervision; higher plot ratios under-represent GEDI and degrade the structural predictions that feed Lphys . A symmetric ratio gives the encoder balanced gradient from both label structures every step, and is the ratio used throughout. Source-balanced sampling does not affect propensity calibration. The source-balanced sampler operates at the patch level, balancing how often GEDI patches and plot patches enter a batch. The propensity π̂ and the AIPW correction operate at the pixel level within a patch, modelling whether a given variable k at pixel (i, j) is observed conditional on the EO covariates x. These two distributions are independent: changing how often patches of each source enter a batch does not change the within-patch distribution of which pixels carry plot labels, which is what the BiasHead is trained to predict. The MNAR mechanism we correct for is therefore the surveyor accessibility bias inside each plot-supervised patch, not the source imbalance between GEDI and plot patches.
D
Implementation and training
The model is implemented in PyTorch [Paszke et al., 2019] and trained with AdamW [Loshchilov and Hutter, 2017], with peak learning rate 5.0e-4 and weight decay 1.0e-4. The learning-rate schedule is cosine annealing with linear warmup over the first 100,000 training steps. We use a batch size of B = 32 (split 50/50 between GEDI and plot samples by the source-balanced sampler; see Appendix C). One epoch corresponds to ⌊NG /(B/2)⌋ ≈ 125, 000 steps. We validate every 500 training steps; the best checkpoint is selected by AGB RMSE on a hold out validation set, with early stopping patience of 250 validation checks (≈ one epoch). Training runs for a maximum of 100 epochs (∼12,500,000 steps); typical convergence is reached within 50 epochs. Full StruMPL training takes approximately 40 hours per site on a single NVIDIA L4 GPU. The architecture, loss components, hyperparameters, and training recipe are identical across the two sites; only the choice of gϕ varies (with vs. without WD). The forward pass conditionally skips the physics, imputation and propensity modules when their loss weights are zero, so ablations run at the cost of a plain multi-task regressor without any code branching. end Loss weight values. We use λstart phys = 0.05, λphys = 0.1, Tphys = 20 epochs, λcons = 0.1, λbias = 0.1, and λimp = 1.0. These were selected on a hold-out validation split via a coarse grid search over each weight independently while keeping the others fixed at their default value.
Numerical stability. The propensity is clamped to [πmin , 1] with πmin = 0.1 to bound the maximum IPW weight at 10×. The physics module is computed in log-space with the structural exponent sp(ŷSD ) · d clamped to [−10, 10] to prevent overflow on early predictions. The physics output is clamped to the ecologically plausible range [0, 2000] Mg/ha before renormalisation. All raw exponent 16
parameters in the physics module are initialised to −4, giving sp(−4) ≈ 0.018 post-softplus and an [ phys in the range 20–25 Mg/ha regardless of input magnitudes. initial AGB Reproducibility. We report headline numbers as mean ± standard deviation over 5 random seeds (42, 123, 456, 789, 1011); ablation rows are single-seed (42). All sources of randomness (PyTorch, NumPy, Python random, and CUDA operations) are seeded in each run.
E
Stop-gradient analysis
We provide complete gradient calculations for the supervised loss Lsup (Equation (4)) under three stop-gradient configurations, demonstrating why both detachments are individually necessary. All 2 gradients are computed pixelwise; the constant factor KHW is absorbed throughout. P 1 Setup. The supervised loss is Lsup = KHW (ŷ − ỹ)2 with pseudo-outcome ỹ = µ + R π̂ (y − µ). At unlabelled pixels (R = 0): ỹ = µ. At labelled pixels (R = 1): ỹ = µ + π̂1 (y − µ). Case 1: both µ and π̂ detached (default). With µ and π̂ treated as constants, the only gradient flow is into ŷ. At labelled pixels: ∂Lsup 1 ∝ 2 ŷ − µ − (y − µ) . ∂ ŷ π̂ At unlabelled pixels: ∂Lsup ∝ 2(ŷ − µ). ∂ ŷ The unlabelled-pixel term provides AIPW’s signal where labels are absent: ŷ is pulled toward the imputation prediction µ, a control-variate signal that does not exist in standard IPW. The imputation m̂ and propensity π̂ are unchanged by Lsup and trained only through their respective losses (Limp , Lbias ). Case 2: only µ detached, π̂ has gradient. Using ∂∂π̂ỹ = − π̂12 (y − µ):
At labelled pixels, the gradient w.r.t. π̂ is non-zero.
∂Lsup 2 ∝ 2 (ŷ − ỹ)(y − µ). ∂ π̂ π̂ Substituting ŷ − ỹ = − π̂1 (y − µ) at the linearization ŷ ≈ µ: ∂Lsup 2 ∝ − 3 (y − µ)2 < 0. ∂ π̂ π̂ This negative gradient on π̂ causes gradient descent to increase π̂, eventually saturating at π̂ → 1. Once π̂ ≈ 1, the IPW correction 1/π̂ has no effect and ỹ = y at labelled pixels, reducing Lsup to naive masked MSE. This is the propensity-collapse failure mode empirically observed in Table 2. Case 3: only π̂ detached, µ has gradient.
ỹ At labelled pixels: ∂∂m̂ = (1 − π̂1 ). So:
∂Lsup ∝ −2(ŷ − ỹ)(1 − π̂1 ) = 2(ŷ − ỹ)( π̂1 − 1). ∂ m̂ ∂L
At unlabelled pixels: ∂ m̂sup ∝ −2(ŷ − µ). Both pull the imputation m̂ toward the regression head ŷ, collapsing m̂ → ŷ. With m̂ ≈ ŷ, the supervised loss at labelled pixels becomes: 2
1 LR=1 sup = (ŷ − ỹ) = − π̂ (y − ŷ)
17
2
= π̂12 (y − ŷ)2 .
This is a 1/π̂ 2 -weighted MSE (strictly worse than standard IPW’s 1/π̂ weighting). Combined with the loss of the AIPW unlabelled-pixel signal (which now contributes zero gradient in expectation as ŷ matches µ), the framework reduces to unstable IPW. This is the imputation-collapse failure mode empirically observed in Table 2. Summary. Both stop-gradients are individually necessary: detaching π̂ prevents the encoder from inflating propensities to escape the MNAR correction; detaching µ preserves the independence between m̂ and ŷ that gives AIPW its control-variate property. Either omission collapses the framework to a degenerate variant of standard MSE or IPW.
F
External baseline evaluation
We compare against two external baselines: Closest published method [Asiyabi et al., 2026]. The [Asiyabi et al., 2026] method was implemented from the published description and trained from scratch on identical training and test splits to those used for StruMPL and the ablations on each site (SNFI for Spain, SEOSAW for Africa). Input channels (Cin = 15), test plot selection, AGB label handling, and evaluation metrics are identical to those used elsewhere in this paper. The training procedure and hyperparameters follow the [Asiyabi et al., 2026] paper. The reported RMSE and bias in Table 1 are computed on the same test pixels as full StruMPL, enabling a fair comparison. This is the closest published method to StruMPL in problem setting and is the more informative of the two external baselines. Off-the-shelf standard reference [Santoro and Cartus, 2024] (ESA CCI Biomass v5.01). The European Space Agency (ESA) Climate Change Initiative (CCI) Biomass map is a global gridded AGB product. Unlike the [Asiyabi et al., 2026] baseline, this is a published map rather than a method we re-ran: we extract the ESA CCI AGB values at the geographic locations of our test plots in both sites and compare them to the corresponding labels from field-plots using the same RMSE and bias metrics as for the trained models. The ESA CCI map is included as a reference point for how an off-the-shelf operational global product compares to site-specific trained models; it is not a like-for-like baseline because it was not trained or calibrated on either of our datasets, and its evaluation is not adjusted for any spatial-resolution or temporal mismatch between the map and the test plots. We report it for transparency and as a benchmark for the order-of-magnitude error a non-customised product delivers, not as a method-level comparison.
G
Physics module variants
The form in equation (5) is one of several differentiable parameterisations we considered. We summarise the alternatives below. All are computed in log-space, all use sp(·) + ε for stability, and all have parameters constrained positive via softplus or exp. Variant power_law allometric (default, this paper) allometric_no_wd mlp
Form
Notes Q
ei
[ = scale · AGB i sp(xi ) α + scale · [H b C c ]SD·d · WDe α + scale · [H b C c ]SD·d [ = MLP([x1 , . . . , xn ]) AGB
n inputs; baseline; n + 1 params 6 params; needs WD 5 params; Spain default expressive but not interpretable
We compare allometric against power_law and mlp as ablations in section 5. allometric_no_wd is structurally allometric without the WD factor and is the default at SNFI; we use the same −4 initialisation strategy across all variants (except for MLP which is initialised randomly upon model definition).
H
Statistical significance of the Spain RMSE improvement
The headline RMSE improvement of full StruMPL over the strongest ablation on Spain (41.9 vs. 43.0 Mg/ha for MTL + Lbias + Limp ) is small relative to the seed-to-seed standard deviation 18
(±0.5 Mg/ha across 5 seeds). To check whether the gap is statistically robust or within sampling noise, we run a paired bootstrap on the Spain test set. Procedure. For each of n = 10,000 bootstrap iterations, we resample the test pixels with replacement (preserving the labelled-pixel distribution and the original sample size) and compute AGB RMSE for both configurations on the resample. The difference ∆ = RMSEStruMPL − RMSEablation is recorded for each iteration. Because both configurations are evaluated on identical resampled pixels, the test is paired: variability in which samples are selected cancels out, and only the model difference contributes to the variance of ∆. Result. The observed RMSE difference is -1.1 Mg/ha (full StruMPL lower). The 95% bootstrap CI on ∆ is [-1.598, -0.522] Mg/ha, entirely below zero. The one-sided p-value (probability of observing ∆ ≥ 0 under the bootstrap distribution) is p < 0.05. We therefore claim the Spain RMSE improvement is statistically significant at the 5% level. The bootstrap shows that full StruMPL’s lower aggregate RMSE on Spain is not an artefact of the test sampling. Why only Spain. The same test is not informative on Africa because the SEOSAW test set contains only 162 plots. At this scale, bootstrap CIs on ∆ are wide enough to admit either sign of the difference for most ablation comparisons, regardless of the underlying truth. We report the Africa point estimates and their 5-seed standard deviations (Table 6) and rely on the larger Spanish test set for the significance claim.
I
Full ablation table
Table 1 and Figure 2 in the main paper report a curated subset of configurations. The complete ablation grid, including single-output baselines, all loss-component combinations, different physic model, naive MSE vs. IPW vs. AIPW, and stop-gradient ablation comparison, is given in Table 5 for all AGB and other structural variables for SNFI and in Table 6 for SEOSAW datasets. Table 5: Full ablation grid on the SNFI (Spain) test set. AGB RMSE in Mg/ha; structural variables in their native units. Configuration
AGB RMSE
AGB bias
H RMSE
C RMSE
SD RMSE
43.7 43.1
0.9 2.1
— —
— —
— —
Joint multi-task variants (no MNAR correction, no physics) MTL, sup only 43.9
1.3
3.1
0.128
459.9
Loss component additions MTL + Lphys MTL + Lbias + Limp MTL + Lcons
0.9 1.5 0.8
3.0 2.9 2.9
0.126 0.124 0.125
459.2 461.4 457.4
3.9 3.1 2.8 ± 0.05
0.125 0.126 0.123 ± 0.001
468.4 469.3 457.9 ± 3.5
Single-output baselines Single-output AGB, naive MSE Single-output AGB + Lbias + Lcons + Limp
43.7 43.0 43.3
Supervised loss form (other components held at full StruMPL) StruMPL with naive masked MSE 43.2 1.7 StruMPL with IPW supervised loss 44.3 0.7 41.9 0.4 StruMPL with AIPW (default) ± 0.5 ± 0.1 Physics module form (other components held at full StruMPL) StruMPL with power_law physics 44.3 StruMPL with mlp physics 44.4
1.0 1.1
3.3 3.0
0.125 0.126
473.2 468.7
Stop-gradient ablation (default = both detached) StruMPL, π̂ not detached 44.2 StruMPL, µ = ŷ not detached 44.6
1.8 1.5
3.9 4.8
0.127 0.123
469.2 469.2
19
Table 6: Full ablation grid on the SEOSAW (Africa) test set. AGB RMSE in Mg/ha; structural variables in their native units. Configuration
AGB RMSE
AGB bias
H RMSE
C RMSE
SD RMSE
WD RMSE
22.5 22.6
5.3 2.9
— —
— —
— —
— —
Joint multi-task variants (no MNAR correction, no physics) MTL, sup only 24.9
1.4
3.9
0.153
97.9
0.111
Loss component additions MTL + Lphys MTL + Lbias + Limp MTL + Lcons
5.4 2.2 1.3
3.9 4.1 4.1
0.152 0.160 0.158
92.6 98.4 108.6
0.112 0.111 0.113
Supervised loss form (other components held at full StruMPL) StruMPL with naive masked MSE 25.3 StruMPL with IPW supervised loss 24.1 20.7 StruMPL with AIPW (default) ± 0.6
1.7 2.6 0.6 ± 0.7
3.9 3.9 3.8 ± 0.1
0.158 0.157 0.153 ± 0.003
94.1 97.8 90.5 ± 3.8
0.115 0.114 0.110 ± 0.007
Physics module form (other components held at full StruMPL) StruMPL with power_law physics 25.5 StruMPL with mlp physics 24.2
7.3 3.4
3.8 3.9
0.151 0.156
112.1 98.6
0.122 0.119
Stop-gradient ablation (default = both detached) StruMPL, π̂ not detached 25.0 StruMPL, µ = ŷ not detached 24.8
3.9 3.6
3.8 3.6
0.155 0.154
92.5 95.7
0.112 0.128
Single-output baselines Single-output AGB, naive MSE Single-output AGB + Lbias + Lcons + Limp
J
24.1 23.6 23.7
Propensity calibration
The AIPW correction relies on π̂ approximating the true propensity. We diagnose the propensity head’s calibration on Spain dataset by binning test pixels by predicted π̂ and computing the empirical labelling rate within each bin (Figure 4).
Figure 4: Propensity calibration on the Spain test set. Predicted π̂ is binned into deciles; the empirical labelling rate within each bin is plotted against the bin midpoint, with the diagonal indicating perfect calibration. Calibration is good at both ends of the range (π̂ ≤ 0.15 and π̂ ≥ 0.85) and weaker in the mid-range, where the model under-predicts the true labelling rate. The monotonic ordering across all bins supports the use of π̂ as a ranking signal in the AIPW correction. Figure 4 shows the predicted propensity π̂ binned into deciles against the empirical labelling rate within each bin, on the Spain test set. Three patterns are worth noting. First, the propensity head is monotonically calibrated: the empirical labelling rate increases with π̂ across all ten bins, from 0.003 at π̂ ≈ 0.05 to 0.991 at π̂ ≈ 0.95. This correct ranking is the property AIPW relies on most directly as the relative ordering of weights determines whether the bias correction re-weights the right 20
pixels. Second, calibration is best at the extremes (π̂ ∈ {0.05, 0.15, 0.85, 0.95}) and weakest in the middle range, where the model under-predicts the true labelling rate (e.g., π̂ ≈ 0.45 vs. empirical rate ≈ 0.88). The practical impact of this mid-range miscalibration is bounded by the propensity clamp πmin = 0.1, which caps the maximum IPW weight at 10× regardless of the true rate; mid-range pixels carry IPW weights between ∼1 and 2, where small calibration errors lead to small weight errors. Third, the BiasHead produces non-trivial spatial variation in π̂ within a single batch: the per-pixel standard deviation of π̂ within a typical test batch is ∼0.13, meaning that the model does not collapse to the degenerate solution in which the BiasHead collapses to a constant π̂ ≈ R̄ and AIPW reduces to a global re-scaling of the loss. Together, these observations support the use of π̂ as the IPW weight in the AIPW correction, not as a maximum-likelihood estimate of the true propensity (which it is not), but as a calibrated ordering of pixels by labelling probability that allows the correction in section 5.2 to focus the supervisory signals on the under-represented regimes of the training data.
K
Cross-source portability
The Africa configuration is structurally different from Spain in three ways: a different biome (dry tropical/savanna vs. Mediterranean), a different KP ({SD, AGB, WD} vs. {SD, AGB}), and a different physics function g. The architecture, optimiser, loss components, weights, and training schedule are identical. Discussion. The Africa results support the portability claim: the unified mask design adapts to the wood density variable availability, the BalancedBatchSampler handles the larger source imbalance (NG /NP ≈ 250 vs. Spain’s ≈ 110) without recipe changes, and the curriculum warmup remains effective with a different physics functional form. We do not claim biome-stratified bias correction at Africa given the test set size (162 standardised 1 ha plots are selected on data-provider recommendation as the most reliable subset of the network).
L
Effect of curriculum warmup
The physics-loss uses a curriculum schedule to increase λphys as the output of the regression head becomes meaningful. Without warmup (λphys at λend phys from step 0), the physics consistency loss penalises the inconsistency of near-random predictions during the first few thousand steps, contributing noise to the gradient. We observe that without this warmup phase, the model convergence is delayed and the final RMSE is higher.
M
When the StruMPL formulation applies
The problem formulation in section 3 is domain-neutral: it refers to inputs x, targets y, an observation mask R, source-determined label subsets Ks , and a known constraint gϕ . The forest-attribute application (section 3, last paragraph) is one example. To clarify the scope of the framework and make the generality claim testable, we list four structural conditions that determine whether StruMPL is the right tool for a given problem. A domain in which all four hold is a candidate. Condition 1: source-determined disjoint partial supervision. Training data come from two or more measurement sources, and the source determines which targets are labelled: Ks ∩ Ks′ = ∅ for s ̸= s′ , with no single sample observing all targets. This rules out standard partial-label settings where any subset of labels may be missing per sample. Condition 2: missingness depends on the latent label. Within at least one source, the labelled subset of locations is not a random sample of the population: P (R = 1 | x, y) ̸= P (R = 1 | x) in some part of the input space. Condition 3: known structural constraint among targets. A function gϕ that links a derived target to other targets is given by domain knowledge, with parameters ϕ either known or learnable. The constraint must be differentiable in the target predictions; for instance: closed-form algebraic 21
relationships, algebraic conservation laws, and parameterised empirical laws. However, pure blackbox constraints (e.g., predictions of a separately trained model) are not what the physics module is designed for. Condition 4: dense or otherwise structured prediction. The target is predicted at many locations per sample (pixels, time steps, mesh nodes), so the propensity field π̂ and the per-location physics residual carry information beyond a scalar. For purely scalar regression, the AIPW correction is still meaningful but the encoder-shared propensity field reduces to a sample-level scalar and the per-location physics constraint becomes a single residual per sample, undermining most of the framework’s value. Candidate application areas. Several settings outside Earth observation plausibly satisfy all four conditions: Hydrology. Streamflow, evapotranspiration, soil moisture, and snow water equivalent are observed by different sensor networks (USGS gauges, FLUXNET towers, SMAP, MODIS, in-situ snow pillows), each covering different catchments and seasons. Snow measurements in particular are MNAR (concentrated where access is feasible). The water balance P = ET + Q + ∆S provides a hard structural constraint among the targets. Battery state estimation. State-of-charge, state-of-health (SoH), capacity, internal resistance, and temperature are observed under different test protocols, with capacity tests typically run only on cells that have not yet failed (an MNAR mechanism on the latent SoH). Coulomb counting and impedance models provide structural relationships among the targets. Multi-modal medical imaging. Different imaging modalities (MRI sequences, CT, PET) are acquired in different patient subsets, with acquisition strongly correlated with disease severity (MNAR). Anatomical-consistency constraints across modalities (though softer than a closed-form law) give a candidate gϕ . However, we note that the constraint condition is weaker than the other candidate cases, since anatomical consistency is empirical rather than physical. Where StruMPL is not the right tool. If labels are randomly missing rather than source-partitioned, a standard partial-label MTL method is sufficient. If observation is MAR rather than MNAR (true for many citizen-science and crowdsourced datasets after geographic stratification), IPW reweighting is unnecessary. If no inter-target structural law is known, the physics module degenerates to a learned regression head with no benefit over additional supervised heads. The framework is targeted at a specific intersection of conditions, not a generic multi-task learner.
22