Bi-FORK: Generative Modeling of High-Dimensional Bifurcating Systems

Anna Zimmel; Fleur Hendriks; Markus Holzleitner; Florian Sestak; Martin Weichselbaumer; Vlado Menkovski; Johannes Brandstetter

cs.LGOct 8, 2026

Bi-FORK: Generative Modeling of High-Dimensional Bifurcating Systems

Authors: Anna Zimmel, Fleur Hendriks, Markus Holzleitner, Florian Sestak, Martin Weichselbaumer, Vlado Menkovski, Johannes Brandstetter

Organizations: ELLIS Unit, LIT AI Lab, Institute for Machine Learning, JKU Linz · Department of Mechanical Engineering, Eindhoven University of Technology · DIFFER – Dutch Institute for Fundamental Energy Research · Department of Mathematics and Computer Science, Eindhoven University of Technology · Mistral, Paris

Abstract

Bifurcations are ubiquitous in physical systems, from structural buckling to fluid and climate dynamics, yet they remain largely unexplored in deep learning. At a symmetry-breaking bifurcation, a single input admits multiple equally valid solutions, violating the one-to-one assumption underlying most learned physical surrogates. We introduce Bi-FORK, a generative framework for learning these one-to-many solution maps in high-dimensional systems. Bi-FORK generates complete trajectories through latent flow matching, preserving space and time coherence, and uses repulsion-guided sampling to recover distinct solution branches in a single amortized pass. We evaluate Bi-FORK on buckling beams, mechanical metamaterials, and Allen-Cahn phase separation, spanning continuous, discrete, and field-valued bifurcations with discretizations up to 260,000 points. Bi-FORK recovers the multimodal solution structure while scaling several orders of magnitude beyond prior approaches, opening generative modeling to high-dimensional bifurcating physical systems.

Figures & tables

Figure 1 : Learning bifurcating solution distributions. Left: Symmetry breaking results in multiple valid trajectories for the same control parameter and initial state. Center: Latent flow matching jointly samples trajectories and resolves distinct branches in one forward pass. Right: Unlike mode-averaging baselines, Bi-FORK recovers the target distribution.

Figure 2 : Exemplary pitchfork bifurcation for the buckling metamaterials use case with increasing control protocol $\bm{\ell}$

Dataset	$G_{\mathbf{c}}$	$H_{\mathbf{c}}$	$\mathcal{M}_{\mathbf{c}}$	$\|\mathcal{M}_{\mathbf{c}}\|$
Beam3D	$\mathrm{O}(2)$	$\mathbb{Z}_{2}$	circle of directions $\phi$	$\infty$
Metamaterials	wallpaper subgroup fixing $\bm{\ell}$	pattern stabilizer	translated patterns	$1,2,4$
Allen-Cahn	domain syms $\times\{u\mapsto\pm u\}$	morphology symmetries	orbits of morphologies	varies

Table 1 : Symmetry structure per dataset: group

G_{\mathbf{c}}

of the condition, isotropy

H_{\mathbf{c}}

of an output state, resulting mode set

\mathcal{M}_{\mathbf{c}}\cong G_{\mathbf{c}}/H_{\mathbf{c}}

(union of orbits for Allen-Cahn).

Figure 3 : Encoder-decoder structure for point cloud datasets. Input tokens (reference positions $\mathbf{X}^{1}$ and displacements) are encoded into latent representation by cross-attending to $L$ learned latent query tokens. The decoder reconstructs the displacements conditioned on the reference positions.

Figure 4 : Particle guidance during sampling. Trajectories evolve identically up to bifurcation point

t_{b}

. From

t_{b}

, the repulsion term acts on the masked frames

t\geq t_{b}

, pushing the

K

trajectories apart, yielding a multimodal distribution.

Figure 5 : Beam3D result comparison. Each panel shows an evaluated method, with an illustrative sample (colored) plotted over the actual trajectory (black). Their endpoints are shown against the full population and the ideal uniform dashed circle. Metrics are aggregated over all Beam3D test conditions using five random seeds.

Figure 6 : Mechanical metamaterials result comparison. Ground truth modes and generated samples are shown for representative conditions with $|\mathcal{M}_{\mathbf{c}}|=2$ and $|\mathcal{M}_{\mathbf{c}}|=4$ buckling modes. Colors indicate the mapped ground-truth mode; missing modes are shown in gray. Affine baseline results in mean prediction. Reported MCon-MAE is averaged over the full test set, while JSD and full coverage are reported for $|\mathcal{M}_{\mathbf{c}}|\in\{2,4\}$ . Metrics reported over five random seeds on the test dataset.

Figure 7 : Allen-Cahn result comparison. Ground truth (black open circles) and predictions (teal and pink circles) are shown across a parameter sweep $\mu\in[-0.1,1.0]$ at fixed $\epsilon=0.1$ . The vertical axis shows the final-frame mean value. We evaluate metrics on the test dataset and average over five random seeds; lower values indicate better performance.

Figure 8 : Visual comparison of ground truth vs predicted samples for $\mu=0.86$ and $\epsilon=0.094$ .

Appendix figures & tables21 assets

Supplementary material from the paper’s appendix.

Appendix

Symbol	Meaning	Space
Data and trajectories
$N$ , $n$	number of spatial locations (nodes); node index	$\mathbb{N}$
$T$ , $t$	number of timesteps; timestep index	$\mathbb{N}$
$D$ , $C$	spatial dimension; channels of the response	$\{2,3\}$ , $\{1,D\}$
$\mathbf{X}^{1}$ , $\mathbf{X}^{t}$	reference geometry; deformed coordinates at timestep $t$	$\mathbb{R}^{N\times D}$
$\mathbf{M}^{t}$	response (signal field) at timestep $t$	$\mathbb{R}^{N\times C}$

Appendix

Table 3 : Overview of used symbols and notations. Superscripts

t

index trajectory frames, subscripts

\tau

flow time, and parenthesized superscripts

(k)

generated samples. Dataset-specific quantities are defined where they are used ( Sections 3 and C ).

Parameter	Range	Add. Information
$N_{\mathrm{seg}}$	3–11 nodes (2–10 segments)	per trajectory
$L_{i}^{\mathrm{beam}}$	log-uniform[ $0.5,2.0$ ]	per segment
$k_{i}^{\mathrm{ax}}$	log-uniform[ $0.5,2.0$ ]	per segment
$k_{i}^{\mathrm{rot}}$	log-uniform[ $0.5,2.0$ ]	per node
$\phi$	uniform on [ $0,2\pi$ )	per trajectory

Appendix

Table 4 : Sampling distributions across used parameters.

Figure 9 : Example trajectories from the 3D beam dataset with different buckling directions $\phi$ for diverse node counts over timestep $t$ .

Split	Trajectories	Share (%)
Training	697	70
Validation	150	15
Test	150	15
Total	997	100

Appendix

Table 5 : Summary of 3D beam dataset splits.

Panel of Figure 10 — Figure 10 : Example mechanical metamaterials trajectories with distinct post-buckling mode structures.

Split	1 mode	2 modes	4 modes	Trajectories	Share (%)
Training	6,529	1,777	262	8,568	70
Validation	1,408	377	51	1,836	15
Test	1,391	354	91	1,836	15
Total	9,328	2,508	404	12,240	100

Appendix

Table 6 : Summary of dataset splits.

Figure 11 : Allen-Cahn example for two distinct $\epsilon,\mu$ settings.

Split	Conditions	Simulated trajectories	Candidate modes	Share (%)
Training	700	2,800	5,600	70
Validation	150	600	1,200	15
Test	150	600	1,200	15
Total	1,000	4,000	8,000	100

Appendix

Table 7 : Summary of the Allen-Cahn dataset splits. Each condition contains four simulated trajectories and up to eight modes after analytic sign augmentation.

Parameter	Beam3D	Metamaterials	Allen-Cahn
Architecture	Perceiver AE	Perceiver AE	ViT3D + spatial token bottleneck
Params	$1.7\times 10^{6}$	$1.7\times 10^{6}$	$6.1\times 10^{6}$
Input shape	node states ( $N\times C$ )	point set (17 groups, 8 node types)	$64{\times}64{\times}64{\times}1$ voxel grid
Patch size	–	–	$8{\times}8{\times}8$ (512 tokens)
Latent shape ( $L\times D_{z}$ )	$6\times 128$	$128\times 128$	$4{\times}4{\times}4$ ( $64\times 384$ )
Encoder blocks (cross/self-attn)	$1\,/\,1$	$1\,/\,1$	ViT depth $6$ (heads 6, embed 192)

Appendix

Table 8 : Hyperparameters of the autoencoder (first stage) used to produce the reported outputs, per dataset.

Parameter	Beam3D	Metamaterials	Allen-Cahn
Params (non-AE)	$7.1\times 10^{6}$	$17.2\times 10^{6}$	$31.8\times 10^{6}$
in_dim / cond_dim	$128\,/\,128$	$128\,/\,128$	$384\,/\,192$
Hidden size	$256$	$512$	$512$
Depth	$5$	$4$	$8$
Heads	$8$	$16$	$16$
MLP ratio	$2.0$	$2.0$	$2.0$

Appendix

Table 9 : Hyperparameters of the latent-space (second-stage) flow model used to produce the reported outputs, per dataset.

Metric	Categories	Tests
Mode Coverage JSD	mapping to $\mathcal{M}_{\mathbf{c}}$ ground truth branches	Are all solution branches reached equally often?
Translation JSD	grid shifts, coarsened to 8 bins per axis ( $8^{3}=512$ )	Given one ground truth sample, how often are periodic shifts recovered?
Rotoreflection JSD	48 octahedral elements	Given one ground truth sample, is the recovered orientation/reflection uniform?

Appendix

Table 10 : Summary of the Allen-Cahn JSD metrics, including their categories and tests.

Table 11 : Beam3D prediction performance. Lower MCon-MAE, JSD, and rejection rate indicate better reconstruction accuracy, better agreement with the target mode distribution, and fewer physically invalid rollouts. (mean

\pm

std over five evaluation seeds)

Table 12 : Result comparison of Bi-FORK (including w/o particle guidance) against the affine baseline, evaluated with MCon-MAE, JSD, and coverage for each mode count. (mean

\pm

std over five evaluation seeds)

Table 13 : Result comparison of Bi-FORK (including w/o particle guidance), evaluated with relaxed MCon-MAE and three JSD metrics. (mean

\pm

std over five evaluation seeds)

Dataset	Accuracy	Precision	Recall	F1	Balanced acc.
Beam3D	$0.989\pm 0.000$	$0.988\pm 0.000$	$1.000\pm 0.000$	$0.994\pm 0.000$	$0.950\pm 0.001$
Metamaterials	$0.970\pm 0.000$	$0.940\pm 0.002$	$0.871\pm 0.002$	$0.904\pm 0.002$	$0.930\pm 0.001$

Appendix

Table 14 : Bifurcation-head classification on the complete test splits. Reporting the mean and standard deviation over five seeds.

Method	Samples	Conditions	Seconds Taken
GeoTDM	180	1	$7766.36$
STFlow	180	1	$72.32$
Bi-FORK w/o PG	180	1	$\mathbf{2.94}$
Bi-FORK	180	1	$3.95$

Appendix

Table 15 : Wall-clock time taken to generate 180 samples of the Beam3D dataset over four random seeds.

Panel of Figure 12 — Figure 12 : Beam bending process: showing predicted (green) vs actual (black) trajectory evolution for 180 samples. Node count ranging from $3$ to $11$ nodes

Figure 13 : Sample comparison for increasing $\mu$ and $\epsilon$ , showcasing the diversity of potential solutions compared to the ground truth.

Panel of Figure 14 — Figure 14 : Visualization of one-mode samples generated by Bi-FORK with $K=16$ . The ground truth is compared with the closest predicted mode.

Panel of Figure 15 — Figure 15 : Visualization of two-mode samples generated by Bi-FORK with $K=16$ . The ground truth is compared with the closest predicted mode.

Figure 16 : Visualization of a four-mode sample generated by Bi-FORK with $K=16$ . The ground truth is compared with the closest predicted mode.

Explore similar work

Branched Optimal Transport Amortization
Sep 14, 2026Semyon Semenov, Viktor Kovalchuk, Meir Roketlishvili +4Flow Matching Generative Modeling
MBZUAI
Bifurcation Models: Learning Set-Valued Solution Maps with Weight-Tied Dynamics
May 8, 2026Caleb Jore, Jialin LiuDynamical Systems
OOD Generalization as a Bifurcation Problem
Sep 27, 2026Nguyen-Thanh-Luong Doan, Quang-Vu Nguyen, Tang-Phu-Quy Le +1Diffusion Model Guidance OOD Generalization
The University of Danang - Vietnam-Korea University of Information and Communication Technology

Sep 14, 2026, cs.LGJ/K move · Enter open · S save

Branched Optimal Transport Amortization

Semyon Semenov, Viktor Kovalchuk, Meir Roketlishvili, Albert Baichorov+3

MBZUAI

Flow Matching Generative Modeling Flow-Based Generative Modeling Optimal Transport

Open paper AskPDF arXiv

Methods of Branched Optimal Transport (BOT) mimic the economy and efficiency of natural tree-like structures, such as those found in rivers and biological systems. These methods are widely applicable for designing efficient networks in society, from river basins and blood vessels to mail and gas distribution systems. However, they remain understudied in the context of designing deep generative models, particularly at a large scale. Standard continuous-time generative models, such as the flow matching approach, fail to capture the inherent hierarchical and branching patterns present in real-world data. Current models provide no mechanism for flows to merge or share pathways to minimize total transport cost. Inspired by the "economy of scale" principle in BOT, we introduce a novel, scalable branched flow-matching algorithm designed to solve the branched optimal transport problem in high dimensions. Our method adapts the Benamou-Brenier continuous-time optimal transport formulation to learn branched generative flows. These flows allow probability mass to aggregate along common pathways before branching out to diverse targets. Parametrized by neural networks, our method effectively learns complex branched generative processes. We demonstrate its effectiveness on challenging high-dimensional tasks in biology and image generation.

Abstract

Explore similar work

Branched Optimal Transport Amortization

Bifurcation Models: Learning Set-Valued Solution Maps with Weight-Tied Dynamics

OOD Generalization as a Bifurcation Problem