@@ -307,6 +307,12 @@ def _compile_partition_edge_block_ids(self) -> None:
         # shapes at `batch_size = 1` gives sizes in per-batch units directly.
         unit_base, ns2unit = 0, dict()
         for ns in self.nodes:
+            # Nodes that share storage share their offsets too, and must not advance the cursor
+            owner = ns.external_params.storage_owner(ns)
+            if owner is not ns and owner in ns2unit:
+                ns2unit[ns] = ns2unit[owner]
+                continue
+
             slots = []
             for shape in ns.external_params.storage_shapes(ns, 1):
                 slots.append(unit_base)
 
@@ -633,16 +633,30 @@ def _external_params_layout(self, batch_size: int):
         if layout is not None:
             return layout
 
+        # Two passes so that a node may share another node's slots (see
+        # `ExternalSumParams.storage_owner`): allocate for the owners first, then alias the rest onto
+        # them. Aliased nodes still appear in the map, so the layers and kernels need no special case --
+        # they look up their own `ns` and simply find the same memory.
         offset = 0
         ns2slots = dict()
         for ns in self.external_params_nodes:
+            if ns.external_params.storage_owner(ns) is not ns:
+                continue
+
             slots = []
             for shape in ns.external_params.storage_shapes(ns, batch_size):
                 slots.append((offset, tuple(shape)))
                 offset += math.prod(shape)
 
             ns2slots[ns] = slots
 
+        for ns in self.external_params_nodes:
+            owner = ns.external_params.storage_owner(ns)
+            if owner is not ns:
+                assert owner in ns2slots, \
+                    f"`storage_owner` of {ns} is not itself an external-parameter node of this circuit."
+                ns2slots[ns] = ns2slots[owner]
+
         layout = (offset, ns2slots)
         self._external_params_layouts[batch_size] = layout
 
@@ -714,9 +728,16 @@ def _stage_external_params(self, kwargs: dict, batch_size: int) -> None:
         torch._foreach_copy_(dsts, srcs)
 
         # Nodes the caller did NOT supply keep whatever the buffer held; drop them so a layer sees
-        # exactly the set that was staged
+        # exactly the set that was staged.
+        #
+        # Membership is by STORAGE GROUP, not by node: when several nodes share one set of tensors (see
+        # `ExternalSumParams.storage_owner`) the caller supplies them once, and every node in that group
+        # is staged. Comparing node-by-node would keep only the one that was named and silently turn the
+        # rest into plain sum layers -- correct-looking output with the correction applied to one
+        # timestep out of many.
+        supplied_groups = {ns.external_params.storage_owner(ns) for ns in ns2tensors}
         for ns in list(views.keys()):
-            if ns not in ns2tensors:
+            if ns.external_params.storage_owner(ns) not in supplied_groups:
                 del views[ns]
 
         kwargs[EXTERNAL_PARAMS_KWARG] = views
 
@@ -42,6 +42,18 @@ def get_signature(self) -> str:
         """
         raise NotImplementedError()
 
+    def storage_owner(self, ns):
+        """
+        The node whose staging slots `ns` uses.
+
+        Identity by default: every node gets its own slots. A parameterization may return a DIFFERENT
+        node to make several nodes share one set of external tensors -- e.g. along a parameter-tying
+        relation, so one factor pair serves every copy of a tied layer. Sharing changes two things for
+        the caller: tensors are supplied once, for the owner, and the gradient returned for the owner is
+        the SUM over every node that shares it.
+        """
+        return ns
+
     def validate_ns(self, ns) -> None:
         """
         Check that `ns` is compatible with this parameterization. Called once, at node construction.
 
@@ -48,7 +48,7 @@ __global__ void lowrank_bw_v_kernel(
         const int64_t* __restrict__ nids, const int64_t* __restrict__ xv,
         float* __restrict__ grad_ext, float* __restrict__ p_lp, float* __restrict__ p_lq,
         int batch_size, int num_eblks, int block_size, int rank,
-        int n_ntiles, int tile_n, long ext_base, int TB) {
+        int n_ntiles, int tile_n, long ext_base, int TB, int accumulate) {
 
     const int re = blockIdx.x;                    // row * num_eblks + edge block
     const int nt = blockIdx.y;                    // node tile
@@ -118,11 +118,21 @@ __global__ void lowrank_bw_v_kernel(
         const long vo = (xv_v + (long)n * rank + r) * batch_size + b;
         const float v = ext[vo];
 
-        // Each (n, r, b) is touched by exactly one thread (nids partition the nodes), so a plain store
-        // is safe and avoids reading a gradient buffer that is as large as `V` itself. Kept in the
-        // exponent rather than forming `f` first, so a tiny flow does not lose the factor.
+        // Kept in the exponent rather than forming `f` first, so a tiny flow does not lose the factor.
         const float lfv = lf + v;
-        grad_ext[vo] = safe_exp_diff(lfv + lw, lt) - safe_exp_diff(lfv + la, lz);
+        const float g = safe_exp_diff(lfv + lw, lt) - safe_exp_diff(lfv + la, lz);
+
+        // A plain STORE is correct when this node owns its factors: `nids` partition the nodes, so each
+        // (n, r, b) is touched exactly once, and it avoids reading a gradient buffer as large as `V`.
+        // When several nodes SHARE one factor pair, every one of them contributes to the same gradient,
+        // so it has to accumulate. Layers run sequentially and each element is still written once per
+        // layer, so a plain read-modify-write suffices -- no atomics -- and the buffer is zeroed once
+        // per backward pass, which is the precondition this relies on.
+        if (accumulate) {
+            grad_ext[vo] += g;
+        } else {
+            grad_ext[vo] = g;
+        }
 
         lse_add(acc_p[w], safe_log_diff(lfv, lt));
         lse_add(acc_q[w], safe_log_diff(lfv, lz));
@@ -181,7 +191,7 @@ __global__ void lowrank_bw_u_kernel(
         const int64_t* __restrict__ cids, const int64_t* __restrict__ xu,
         float* __restrict__ grad_ext, float* __restrict__ element_flows,
         int batch_size, int num_edges, int num_eblks, int ch_block_size, int rank,
-        int tile_c, long ext_base, int TB) {
+        int tile_c, long ext_base, int TB, int accumulate) {
 
     const int re = blockIdx.x;
     const int ct = blockIdx.y;
@@ -214,8 +224,16 @@ __global__ void lowrank_bw_u_kernel(
         const float u = ext[uo];
 
         // The positive term is both dLL/dU's first part and this child's share of the flow
-        grad_ext[uo] = expf(u + e + lp) - expf(u + lq);
-        lse_add(log_corr, u + e + lp);
+        const float log_t1 = u + e + lp;                   // LOG space -- `log_corr` accumulates logs
+        const float g = expf(log_t1) - expf(u + lq);
+
+        if (accumulate) {                          // see the note in pass A
+            grad_ext[uo] += g;
+        } else {
+            grad_ext[uo] = g;
+        }
+
+        lse_add(log_corr, log_t1);
     }
 
     // Atomic because several parent blocks may share a child block: unlike `node_mars` in the forward,
@@ -266,7 +284,7 @@ void lowrank_backward(torch::Tensor node_flows, torch::Tensor element_flows,
                       torch::Tensor p_lp, torch::Tensor p_lq,
                       torch::Tensor log_p, torch::Tensor log_q,
                       int64_t block_size, int64_t ch_block_size, int64_t rank, int64_t ext_base,
-                      int64_t tile_n, int64_t tile_c, int64_t tb) {
+                      int64_t tile_n, int64_t tile_c, int64_t tb, bool accumulate) {
 
     const int batch_size = node_flows.size(1);
     const int rows = xu.size(0);
@@ -285,7 +303,7 @@ void lowrank_backward(torch::Tensor node_flows, torch::Tensor element_flows,
         nids.data_ptr<int64_t>(), xv.data_ptr<int64_t>(),
         grad_ext.data_ptr<float>(), p_lp.data_ptr<float>(), p_lq.data_ptr<float>(),
         batch_size, num_eblks, (int)block_size, (int)rank,
-        n_ntiles, (int)tile_n, (long)ext_base, (int)tb);
+        n_ntiles, (int)tile_n, (long)ext_base, (int)tb, accumulate ? 1 : 0);
 
     const int n_outputs = rows * num_eblks * (int)rank * batch_size;
     lowrank_bw_pq_reduce_kernel<<<bw_cdiv(n_outputs * 32, 256), 256, 0, stream>>>(
@@ -301,5 +319,5 @@ void lowrank_backward(torch::Tensor node_flows, torch::Tensor element_flows,
         cids.data_ptr<int64_t>(), xu.data_ptr<int64_t>(),
         grad_ext.data_ptr<float>(), element_flows.data_ptr<float>(),
         batch_size, num_edges, num_eblks, (int)ch_block_size, (int)rank,
-        (int)tile_c, (long)ext_base, (int)tb);
+        (int)tile_c, (long)ext_base, (int)tb, accumulate ? 1 : 0);
 }
@@ -288,7 +288,7 @@ void lowrank_backward(torch::Tensor node_flows, torch::Tensor element_flows,
                       torch::Tensor p_lp, torch::Tensor p_lq,
                       torch::Tensor log_p, torch::Tensor log_q,
                       int64_t block_size, int64_t ch_block_size, int64_t rank, int64_t ext_base,
-                      int64_t tile_n, int64_t tile_c, int64_t tb);
+                      int64_t tile_n, int64_t tile_c, int64_t tb, bool accumulate);
 
 void lowrank_shift_logz(torch::Tensor node_mars, torch::Tensor nids, torch::Tensor log_z,
                         int64_t block_size, double sign);