Reverse-mode autodiff: order-2 tape slows order-1 gradients, rewound tapes keep allocating blocks, rvar copies silently share a node

Start with the three standalone repros in the issue, built against include/boost/math/differentiation/autodiff_reverse.hpp and reverse_mode_autodiff_memory_management.hpp. Trace the order-2 tape's flat_linear_allocator construction and clear() path to see where placeholder rvar leaves get created, then the emplace_back and emplace_back_n block-allocation logic for the rewind leak. Done means the leaf count stays at zero after the order-2 tape is created, RSS stays flat across rewind cycles, and the copy-aliasing behavior is either documented in doc/differentiation/autodiff_reverse.qbk or changed, with the maintainers choosing which.

I hit three problems in autodiff_reverse.hpp while using it in an optimizer. Each one is a separate standalone program compiled with clang++ -std=c++17 -O2 -DBOOST_MATH_STANDALONE on Apple M-series (arm64), at develop 6dc59aa15. All three start with:

#include <boost/math/differentiation/autodiff_reverse.hpp>
#include <iostream>
#include <vector>
namespace rd = boost::math::differentiation::reverse_mode;
using rv = rd::rvar<double, 1>;

1. Touching the order-2 tape puts 262144 leaves on the order-1 tape

Once get_active_tape<double, 2>() has been called, the order-1 tape holds 262144 extra leaf nodes. While they stay on the tape, every order-1 grad() costs about 0.7 ms more, even for a function of two variables. Each time below is one loop iteration: build the graph, call grad(), then rewind_to_last_checkpoint().

int main() {
  auto& t1 = rd::get_active_tape<double, 1>();
  (void)rd::get_active_tape<double, 2>();
  std::cout << (t1.end() - t1.begin()) << "\n";   // 262144
}

Cause: the order-2 tape stores its adjoints and derivatives as flat_linear_allocator<rvar<double,1>, 65536> (autodiff_reverse.hpp:83-84). Every block is a std::array of rvar<double,1> built by make_unique (reverse_mode_autodiff_memory_management.hpp:246), so each block runs the rvar() constructor 65536 times and each call adds a leaf to the order-1 tape (autodiff_reverse.hpp:434-437). The allocator constructor allocates one block (memory_management.hpp:250). Then gradient_tape() calls clear() (autodiff_reverse.hpp:125, 132), which frees it and allocates a second one (memory_management.hpp:274-277). That is 2 allocators x 2 blocks x 65536 = 262144 leaves, and each later block of those two allocators adds 65536 more. Every order-1 grad() makes three linear passes, and all three visit the placeholder leaves. zero_grad() runs over all adjoints (autodiff_reverse.hpp:664). Then backward() makes the other two: tape.find(node_) (autodiff_reverse.hpp:579, a std::find_if from begin() at memory_management.hpp:442) and the walk back to tape.begin() (autodiff_reverse.hpp:581).

2. Rewinding a tape doesn't reuse its blocks

A loop that builds a graph larger than one block, then calls rewind_to_last_checkpoint() (or rewind()), uses 8 MiB more memory every iteration, while the node count stays the same:

#include <mach/mach.h>   // resident size on macOS
static double rss_mb() {
  mach_task_basic_info_data_t i; mach_msg_type_number_t c = MACH_TASK_BASIC_INFO_COUNT;
  task_info(mach_task_self(), MACH_TASK_BASIC_INFO, (task_info_t)&i, &c);
  return i.resident_size / 1048576.0;   // MiB
}
int main() {
  auto& tape = rd::get_active_tape<double, 1>();
  rv x(1.0);
  tape.add_checkpoint();
  for (int cycle = 0; cycle <= 200; ++cycle) {
    rv f = x;
    for (int i = 0; i < 100000; ++i) f = f * 1.0000001;  // crosses one 65536 boundary
    f.backward();
    if (cycle % 50 == 0)
      std::cout << cycle << ": nodes " << (tape.end() - tape.begin()) << ", RSS " << rss_mb() << " MiB\n";
    tape.zero_grad();
    tape.rewind_to_last_checkpoint();
  }
}

0: nodes 100001, RSS 17.7 MiB
50: nodes 100001, RSS 418.8 MiB
100: nodes 100001, RSS 819.8 MiB
200: nodes 100001, RSS 1622 MiB

Cause: rewinding only resets total_size_ (memory_management.hpp:282, 300). The next emplace_back that reaches a block boundary calls allocate_buffer() (memory_management.hpp:319-320, 336-337) even if data_[buffer_id()] already exists. The new block goes onto the end of data_ and is never used, because indexing still goes through data_[i / buffer_size]. That is one wasted block per allocator per boundary crossing, and the tape has four allocators. One block of each is 65536 x (8 + 8 + 104 + 8) bytes = 8 MiB, where 104 bytes is sizeof(gradient_node<double,1>). emplace_back_n has the same problem (memory_management.hpp:367, 400). It also constructs the elements for the next block in data_.back() (line 368, 401), which is not the block the returned iterator reads. A fix is for emplace_back to allocate only when buffer_id() == data_.size(), and for emplace_back_n to allocate only when block bid + 1 doesn't exist yet and to construct its elements there. emplace_back_n also takes the crossing branch when a run ends exactly at a block boundary (iid + n < buffer_size, lines 353 and 386), which allocates an extra block. And it can only cross one boundary, so a runtime n larger than a block would need a loop.

3. Copies of an rvar share its node, and the docs don't say so

The copy constructor and copy assignment are defaulted (autodiff_reverse.hpp:453-454), so a copy refers to the same tape node and adjoint. That may be the intended handle semantics, but doc/differentiation/autodiff_reverse.qbk only lists the copy and assignment signatures. It never says that copies alias. The behavior is hard to spot because operator=(RealType) only makes a new leaf when node_ == nullptr (autodiff_reverse.hpp:445-450). Copies that are later given different values still share one node, so they get one combined gradient:

int main() {
  std::vector<rv> v(3, rv(0.0));
  v[0] = 1.0; v[1] = 2.0; v[2] = 3.0;
  rv f = v[0] * v[0] + 10 * v[1] + 100 * v[2];
  auto g = rd::grad(f, &v[0], &v[1], &v[2]);
  std::cout << g[0] << " " << g[1] << " " << g[2] << "\n";  // 112 112 112, expected 2 10 100
}

With std::vector<rv> v(3); instead, which default-constructs separate leaves, the result is the correct 2 10 100. The same happens with rv a = 1.0; rv b = a; b = 5.0;: grad(a + 2*b, &a, &b) gives (3, 3) instead of (1, 2). Under shared-node semantics these numbers are consistent, so this may be a documentation issue rather than wrong differentiation. Writing v[i] = rv(value) already gives independent leaves. Making operator=(RealType) always create a new leaf would change the contract: a leaf created after a checkpoint would be dropped by rewind_to_last_checkpoint(). So the smaller change may be to document that copies alias, with the vector(n, value) pitfall as the example. The library optimizers update parameters through get_value(), so either choice leaves them working.