perf_aad.cpp

// SPDX-License-Identifier: Apache-2.0
// The real AAD value proposition: the entire gradient (every sensitivity) from
// one reverse sweep, at a cost independent of the number of inputs. We price a
// book of N vanilla calls on N underlyings and compute all N deltas.
//
//   AAD  : one forward evaluation + one reverse sweep      -> O(N)
//   Bump : 2N full portfolio revaluations (central diff)   -> O(N^2)
//
// For a small book the tape overhead makes bump faster; as N grows AAD wins
// decisively. Both are printed so the crossover is visible and honest.

#include <chrono>
#include <cmath>
#include <cstdio>
#include <vector>

#include "radonlab/aad.hpp"
#include "radonlab/analytics.hpp"

using clk = std::chrono::steady_clock;
using namespace radonlab;

namespace {
constexpr double R = 0.03, SIG = 0.25, T = 1.0, K = 100.0;

double portfolio_plain(const std::vector<double>& spots) {
    double v = 0.0;
    for (double s : spots) v += analytics::vanilla_call(s, K, R, SIG, T, 0.0, "price");
    return v;
}
}  // namespace

int main() {
    std::printf("All deltas of an N-option book: AAD (1 sweep) vs bump (2N revaluations)\n\n");
    std::printf("%6s %14s %14s %10s %14s\n", "N", "AAD (s)", "bump (s)", "speedup", "max|grad diff|");

    std::vector<double> adj;
    for (std::size_t N : {1u, 4u, 16u, 64u, 256u, 1024u}) {
        std::vector<double> spots(N);
        for (std::size_t i = 0; i < N; ++i) spots[i] = 90.0 + 20.0 * static_cast<double>(i) / N;

        const std::size_t reps = std::max<std::size_t>(1, 2'000'000 / (N + 1));

        // ---- AAD gradient ----
        std::vector<double> grad_aad(N);
        auto t0 = clk::now();
        volatile double sink = 0.0;
        for (std::size_t rep = 0; rep < reps; ++rep) {
            ad::tape().clear();
            std::vector<ad::advar> S;
            S.reserve(N);
            for (double s : spots) S.emplace_back(s);  // record spot leaves
            ad::advar value(0.0);
            for (const ad::advar& si : S) value = value + bs_call_advar(si, K, R, SIG, T);
            ad::tape().adjoints_into(value.index(), adj);
            for (std::size_t i = 0; i < N; ++i)
                grad_aad[i] = adj[static_cast<std::size_t>(S[i].index())];
            sink += grad_aad[0];
        }
        const double aad_s = std::chrono::duration<double>(clk::now() - t0).count();

        // ---- Bump gradient (central differences, full revaluation) ----
        std::vector<double> grad_bump(N);
        t0 = clk::now();
        for (std::size_t rep = 0; rep < reps; ++rep) {
            for (std::size_t i = 0; i < N; ++i) {
                const double h = 1e-4 * spots[i];
                std::vector<double> up = spots, dn = spots;
                up[i] += h;
                dn[i] -= h;
                grad_bump[i] = (portfolio_plain(up) - portfolio_plain(dn)) / (2 * h);
            }
            sink += grad_bump[0];
        }
        const double bump_s = std::chrono::duration<double>(clk::now() - t0).count();
        (void)sink;

        double maxdiff = 0.0;
        for (std::size_t i = 0; i < N; ++i)
            maxdiff = std::max(maxdiff, std::abs(grad_aad[i] - grad_bump[i]));

        std::printf("%6zu %14.4f %14.4f %10.2fx %14.2e\n", N, aad_s, bump_s, bump_s / aad_s, maxdiff);
    }

    std::printf(
        "\nAAD returns the whole delta vector from a single reverse sweep, so its\n"
        "cost grows with the book size once, not once per sensitivity. This is the\n"
        "regime that matters for a risk book with many underlyings.\n");
    return 0;
}