Scalable Post-Nonlinear Causal Discovery via Likelihood-Based Order Learning

Post-nonlinear (PNL) structural causal models allow for flexible nonlinear mechanisms but remain difficult to learn in multivariate and high-dimensional settings. The few existing multivariate PNL methods primarily estimate a causal order through recursive sink identification and scale poorly. We present SCORPIO, a method for learning multivariate PNL DAGs under causal sufficiency. SCORPIO parameterizes each inverse outer transformation as a strictly monotone sigmoid normalizing flow, yielding a tractable node-wise likelihood that supports several noise distribution families. We estimate a topological order by greedily adding edges based on local likelihood improvements, which requires updating only the affected node model, avoiding retraining of the full system model. Candidate edges are subsequently pruned via marginal independence tests between scalar residuals. Across synthetic experiments varying the noise distribution, parental-mechanism complexity, graph size, and graph density, SCORPIO achieves the lowest order divergence and structural Hamming distance in most settings. Under a 24-hour time limit, SCORPIO is the only evaluated multivariate PNL method that completed experiments beyond $100$ variables, with results reported for up to $500$ variables. It also achieves the lowest order divergence and competitive DAG-recovery performance on a semi-synthetic manufacturing benchmark and the Sachs protein dataset.