Finding Conservation Laws of Large Dynamical Systems with Tasks and Futures:
A Case Study in Utilizing Dynamic Data Dependencies

We adopt the definition and implementation of futures described in [5] and [6]. Futures act as placeholders for values, featuring distinct writer and reader sides. Specifically, each placeholder is associated with a single promise (used to write the value) and one or more futures (used to read the value). The promise is fulfilled exactly once; upon fulfillment, all associated futures become ready immediately.

Promises and futures are identified by their name and type, and can be stored in standard data structures such as arrays. The element type is unrestricted, enabling advanced constructs such as futures of futures. When passed as parameters, futures are duplicated as needed, while promises remain unique and are always moved (and never copied). Moving a promise transfers the responsibility for filling it to the recipient task. A task is ready to execute only when all futures it depends on have been fulfilled.

Regarding memory management, our implementation utilizes a reference-counting scheme. When all futures referencing a specific placeholder go out of scope, the placeholder and its underlying memory are automatically deallocated, providing a limited form of automatic memory reclamation.

We will use the general notion that spawning a task will return a future for the result of the task. A task can also be spawned into a promise, thereby binding the result of that task to a pre-existing promise. The code in Listing 1 demonstrates how recursive computations and dependencies are expressed using the promise/future model.

• add:

  When used to create a task, takes two integers, computing their sum, and fulfills the output promise p via p.set(i+j).

• fib:

  Computes the $n$-th Fibonacci number $F(n)$ and fulfills promise p:

  • –

    Base case: For $n<2$, fulfills p immediately with $1$.

  • –

    Recursive case:

    1. 1.

       Spawn fib(n-1), obtaining future f1 for $F(n-1)$.

    2. 2.

       Create promise p2 for $F(n-2)$ and obtain its future f2.

    3. 3.

       Spawn add with f1, f2, and output promise p (scheduled only after both futures are ready)

    4. 4.

       Compute $F(n-2)$ directly via fib(n-2, p2) (synchronous in this example).

The application of futures to $\mathcal{H}$-matrix arithmetic, as detailed in [5], leverages the recursive nature of the data structure. Non-elementary tasks recursively spawn other non-elementary tasks as they traverse the matrix hierarchy, eventually reaching the leaf level where elementary tasks are executed. Elementary tasks perform the actual, numerical work; in our implementation, this will be done by an external library. Futures are used to track dependencies throughout this process.

As tasks execute, several pieces of information must be communicated between them to manage dependencies and data flow:

• (i1)

  Has the task been completed? For non-elementary tasks, this means that all child tasks have been spawned, whereas for elementary tasks it indicates that the associated matrix block has been updated. The cases are combined, since a (non-elementary) task generally does not know whether another task is elementary or not.

• (i2)

  Is the update of the associated matrix block complete? (This is only relevant for non-elementary tasks, where (i2) differs from (i1))

• (i3)

  Contents of the (updated) matrix block (This is only relevant for elementary tasks)

• (i4)

  Futures of all child tasks (This is only relevant for non-elementary tasks)

To facilitate this communication, we define the recursive node structure in Listing 2. The rec_matrix_node struct contains a data future for storing the memory address of leaf blocks and a boolean flag leaf to distinguish between internal and leaf nodes. The matrix_t type represents the underlying matrix type used by the external library mentioned above. The array of futures subnodes contains the futures for child tasks in non-leaf nodes. Finally, the done future signals the completion of the operation associated with the node.

It is important to distinguish the roles of these futures: data becoming ready signifies that the memory address for the block is known (essential for non-elementary tasks to schedule children), whereas done becoming ready signifies that the computation on that block (and all subnodes) has finished. Which task is responsible for setting them depends on the algorithm being implemented.

The primary motivation for our matrix inversion algorithm stems from the problem of identifying conservation laws in dynamical systems.

Consider a dynamical system defined by $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$, where $\mathbf{x}\in\mathbb{R}^{n}$. A conservation law is a scalar function $C(\mathbf{x})$ that remains constant along the trajectories of the system, satisfying $\nabla C(\mathbf{x})\cdot\mathbf{f}(\mathbf{x})=0$.

As proposed by Mebrartie et al. [4], this problem can be formulated as a regression task using Kernel Ridge Regression (KRR). Given a set of trajectory data $\mathcal{D}=\{(\mathbf{x}_{(i,j)})_{i=1}^{N}\}_{j=1}^{M}$, where $N$ is the number of samples per trajectory and $M$ is the number of trajectories, the objective is to find a function $C(\mathbf{z})$ that minimizes the deviation from constancy along each trajectory.

The conservation law $C(\mathbf{z})=\sum_{i=1}^{NM}\alpha_{i}\varphi(\mathbf{z},\mathbf{x}_{i})$ is then expressed as a linear combination of kernel functions $\varphi(\mathbf{x},\mathbf{y})=(\langle\mathbf{x},\mathbf{y}\rangle-c)^{d}$, where $\boldsymbol{\alpha}\in\mathbb{R}^{NM}$ are the coefficients to be determined.

To find the coefficients $\boldsymbol{\alpha}$, one must solve a quadratic optimization problem. This involves the inversion of the regularized kernel matrix $\mathbf{K}-\lambda\mathbf{I}\in\mathbb{R}^{NM\times NM}$, where $\lambda>0$ is a regularization parameter and the symmetric matrix $\mathbf{K}$ is computed from the data points $\mathcal{D}$.

The authors note that execution time is dominated by the repeated inversion of kernel matrices. Since matrix size depends on the number of training samples ($N\cdot M$), these matrices can become very large. Furthermore, identifying optimal parameters via cross-validation requires solving this inverse problem multiple times for different matrix configurations.

The await_delete construct takes an existing future $f$ as input and returns a new future $f^{\prime}$. The promise associated with $f^{\prime}$ is implicit and tied to the lifetime of the original future $f$. Specifically, $f^{\prime}$ is fulfilled only when the reference count of $f$ drops to zero, indicating that all consumers of the original value have released their handles.

From the perspective of a task waiting on $f^{\prime}$, the behavior is identical to waiting on a standard future. The fulfillment of $f^{\prime}$ then signals not just that data is available, but that the exclusive ownership of the underlying memory has been reclaimed. This allows the consumer of $f^{\prime}$ to safely treat the memory as mutable and reuse it for subsequent operations without allocating new buffers.

It is important to note that await_delete can only be invoked once per future, as it assumes that the reference count is sufficient to track exclusivity. Provided this invariant holds, scheduling tasks according to their future-based dependencies (as outlined in Subsection 2.1) guarantees deterministic execution order. However, as with standard parallel floating-point arithmetic, bit-wise reproducibility is not guaranteed due to the non-associative nature of floating-point operations.

Furthermore, the await_delete construct composes safely within arbitrary dependency graphs, including nested structures. Whether $f$ is the result of a simple computation or a complex chain of nested futures, $f^{\prime}$ is fulfilled only when the specific future $f$ becomes unreachable. This local property ensures that await_delete can be applied to any node in a graph without introducing race conditions or deadlocks, provided the write-once constraint of the promise is maintained and await_delete is invoked at most once per future.

We apply the await_delete construct to the problem of inverting a symmetric matrix using a standard block-wise recursive approach. Given a symmetric matrix $A$ partitioned into four blocks, the inverse $A^{-1}$ can be computed using the following block decomposition:

where the intermediate steps are defined as:

Note that each operation overwrites the respective block; the notation above distinguishes the intermediate values ($T$, $S$) to clarify the data flow.

Translating this into the future-model is then straightforward; each overwrite happens on an await_deleted future, where the result is passed to the subsequent operations. Algorithm 1 details this process.

The function Invert takes a recursive node representing a matrix block. If the node is a leaf, we simply invoke an external routine to invert the dense block in-place.

For non-leaf nodes, the algorithm follows the dependencies defined in the block decomposition. We spawn tasks to compute the intermediate components ($T_{11},T_{12},T_{22}$) and tasks for the final inverse blocks ($S_{11},\dots,S_{22}$) sequentially.

Crucially, each step reuses the memory of a specific input block:

• •

  Line 5: We recursively invert the top-left block $A_{11}$. The result $T_{11}$ is stored directly in the memory previously occupied by $A_{11}$.

• •

  Line 6: We compute $T_{12}=A_{21}T_{11}$. The result will overwrite the top-right block $A_{12}$.

• •

  Line 7: We compute the Schur complement $T_{22}=A_{22}-T_{12}A_{21}$. This is a fused multiply-add operation that will overwrite the bottom-right block $A_{22}$.

• •

  Line 8: We recursively invert the Schur complement $T_{22}$ to obtain $S_{22}$, reusing the same memory slot.

• •

  Line 9: We compute $S_{21}=-S_{22}T_{12}$. This will overwrite the bottom-left block $A_{21}$.

• •

  Line 10: We compute $S_{12}$ as the transpose of $S_{21}$. This reuses the memory currently holding $T_{12}$ (which was previously $A_{12}$).

• •

  Line 11: Finally, we compute $S_{11}=T_{11}-S_{21}^{T}S_{12}$. This will overwrite the memory holding $T_{11}$ (which was originally $A_{11}$).

While standard block-wise algorithms overwrite input blocks to save memory, traditional futures would necessitate allocating new memory for every intermediate result, incurring significant overhead. The await_delete construct resolves this by enabling precise memory reuse sequencing. It signals when a block is no longer needed by dependent tasks, allowing its memory to be safely reclaimed for subsequent operations. This creates a pipeline that recycles memory exactly between computational phases, maintaining parallelism while minimizing footprint as described in Subsection 3.1.2.

Note that this algorithm does not make any assumptions about the storage layout of the leaf blocks. Indeed, the dependency management with futures allows the computation to flow downward, ultimately overwriting each (leaf) block without the tasks higher up in the recursion hierarchy needing to know where they are actually stored. By the same token, it also generalizes to both operations on non-square matrices, as well as the inversion of non-symmetric matrices. However, while inverting non-symmetric matrices will require additional working space, utilizing the await_delete construct allows us to bound this overhead. Specifically, we can reuse the memory of the input matrix to store the result (as above), but also of the additional working space, thereby limiting the total auxiliary memory requirement to at most half the size of the matrix.

Figure 1 presents the execution time measurements. Each subfigure corresponds to one matrix dimension, with the three curves representing the different block size ratios. The results confirm that parallelism is effectively exploited, as execution time decreases with increasing core count for all configurations. However, the impact of task granularity is significant. In Subfigure 1A ($2^{14}\times 2^{14}$), we observe that finer granularity (smaller block sizes, corresponding to higher ratios) worsens performance. Specifically, the configuration with a ratio of $512$ (smallest blocks) exhibits notably higher runtimes compared to the coarser configurations. For the larger matrix sizes ($2^{15}$ and $2^{16}$), the performance gap between block sizes $64$ and $128$ (ratios $256$ and $128$) is less pronounced, though the trend of coarser granularity yielding better performance persists.

To better quantify the efficiency of the parallelization, Figure 2 depicts the computed speedup curves. Each subfigure corresponds to a ratio between matrix size and block size. For the largest matrix size ($2^{16}\times 2^{16}$), the speedup curves are nearly identical across all block size ratios (Subfigure 2C). This suggests that for sufficiently large problems, the overhead associated with fine-grained task management becomes negligible compared to the computational work. For the intermediate matrix size ($2^{15}\times 2^{15}$), the optimal performance is achieved at a ratio of $256$ (Subfigure 2B). Here, the finer granularity of ratio $512$ leads to a noticeable degradation in speedup, while the coarser ratio of $128$ performs similarly to the optimal case. Interestingly, for the smallest matrix size ($2^{14}\times 2^{14}$), the speedup at ratio $128$ converges with the performance of the $2^{15}$ matrix at the same ratio, indicating that the task granularity is sufficiently coarse to amortize the overhead even for smaller problem sizes.

The performance degradation at fine granularities points to a runtime bottleneck rather than an algorithmic flaw. We attribute this to the current future implementation, which uses a locking scheme (atomic compare-and-swap) for reference counting. With many small tasks, congestion arises as numerous threads simultaneously attempt to read or fulfill futures, causing significant synchronization overhead. This effect is most pronounced when the task creation-to-computation time ratio is high. As matrix size increases, the computation-to-communication ratio improves, diminishing the relative cost of lock contention and closing the performance gap for the $2^{16}$ case.

To validate our hypothesis that lock contention is the primary bottleneck, we profiled the execution to separate the time spent in actual computation from the time spent managing futures. Figure 3 provides a detailed breakdown of the runtime for matrix sizes $2^{14}$ (Subfigure 3A) and $2^{15}$ (Subfigure 3B), categorizing the time into three components: compute time (actual numerical computation), future management (every operation that touches the shared state of futures), and RTS overhead.

The data clearly illustrates the trade-off between granularity and overhead. For configurations with fine granularity (high ratio $r$), the time spent on future management (red) constitutes a significant portion of the total runtime. This confirms that the cost of managing dependencies dominates when tasks are small. Conversely, as the block size increases (lower ratio $r$), the proportion of time spent on future management drops and the application becomes compute-bound.

Comparing this approach to a traditional threading model like OpenMP would likely show lower overheads for the dependency management. However, such models lack the dynamic flexibility required for our target application, specifically the handling of irregular, runtime-generated dependencies inherent to $\mathcal{H}$-matrix arithmetic and the conservation law problem (Subsections 2.2 and 2.3). Thus, the observed overhead represents the cost of this increased expressiveness.