Performance hit using AAD

[fredrik-fyring]

First, thanks for the outstanding work on enabling AAD in QuantLib.

I have a problem with performance degradation when pricing many swaps with AD enabled. Pricing 500 swaps is taking 3x using AD than without. Pricing 10,000 swaps, however, takes 18x time (36x using Release build). The performance hit scales linearly so the XAD code seems to have O(N^2) time complexity. In other words, each NPV calculation seem to be dependent in some way of the previous.

Is this expected?

Portfolio Size	Time without XAD [ms]	Time with XAD [ms]	Performance Diff
2000	303	2590	9
4000	607	11671	19
6000	899	21119	23
8000	1239	36602	30
10000	1513	55938	37

My example code seems quite straight forward, following the AdjointSwapXAD in QuantLib-Risks. First I build a curve and then I use this curve to price and get sensitivities per trade.

I built the example below with "windows-xad clang-debug" and "windows-xad-clang-release" configuration on Windows 11.

#include <ql/qldefines.hpp>
#if !defined(BOOST_ALL_NO_LIB) && defined(BOOST_MSVC)
#    include <ql/auto_link.hpp>
#endif
#include <ql/indexes/ibor/euribor.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
#include <ql/pricingengines/swap/discountingswapengine.hpp>
#include <ql/quotes/simplequote.hpp>
#include <ql/termstructures/iterativebootstrap.hpp>
#include <ql/termstructures/yield/bootstraptraits.hpp>
#include <ql/termstructures/yield/piecewiseyieldcurve.hpp>
#include <ql/termstructures/yield/ratehelpers.hpp>
#include <ql/termstructures/yieldtermstructure.hpp>
#include <ql/time/calendars/target.hpp>
#include <ql/time/daycounters/actual360.hpp>
#include <ql/time/daycounters/thirty360.hpp>
#include <XAD/XAD.hpp>
#include <chrono>
#include <vector>
#include <iostream>
#include <map>
#include <tuple>
#include <fstream>
#include <filesystem>

using namespace QuantLib;

const int Ndepos = 10;
const int Nfra = 5;

// prepares a set of market quotes for deposits, FRAs, and Swap rates, used for curve bootstrapping
void prepareQuotes(Size maximumMaturity, std::vector<double>& marketQuotes) {
    // setting up market quotes
    // deposit quotes on, tn, sn, sw, 1m, ... , 6m
    for (Size i = 0; i < Ndepos; ++i)
        marketQuotes.push_back(0.0010 + i * 0.0002);
    // fra quotes 1-7, ... , 5-11
    for (Size i = 0; i < Nfra; ++i)
        marketQuotes.push_back(0.0030 + i * 0.0005);
    // swap quotes 1y, ... , maximum maturity
    for (Size i = 0; i < maximumMaturity; ++i)
        marketQuotes.push_back(0.0060 + i * 0.0001);
}

// bootstraps the curve used for the swaps
Handle<YieldTermStructure>
bootstrapCurve(Date referenceDate, const std::vector<Real>& marketQuotes, Size maximumMaturity) {
    // build quotes
    std::vector<ext::shared_ptr<SimpleQuote>> quotes;
    std::transform(marketQuotes.begin(), marketQuotes.end(), std::back_inserter(quotes),
        [](const Real& quote) { return ext::make_shared<SimpleQuote>(quote); });

    std::vector<RelinkableHandle<Quote>> quoteHandles;
    std::transform(quotes.begin(), quotes.end(), std::back_inserter(quoteHandles),
        [](ext::shared_ptr<SimpleQuote>& q) { return RelinkableHandle<Quote>(q); });

    // rate helpers
    std::vector<ext::shared_ptr<RateHelper>> instruments;

    // deposits
    for (Size i = 0; i < Ndepos; ++i) {
        Period matTmp;
        Size fixingDays = 2;
        switch (i) {
        case 0:
            matTmp = 1 * Days;
            fixingDays = 0;
            break;
        case 1:
            matTmp = 1 * Days;
            fixingDays = 1;
            break;
        case 2:
            matTmp = 1 * Days;
            break;
        case 3:
            matTmp = 1 * Weeks;
            break;
        default:
            matTmp = (i - 3) * Months;
            break;
        }
        instruments.push_back(ext::make_shared<DepositRateHelper>(
            quoteHandles[i], matTmp, fixingDays, TARGET(), ModifiedFollowing, false, Actual360()));
    }

    // FRAs
    for (Size i = 0; i < Nfra; ++i) {
        instruments.push_back(
            ext::make_shared<FraRateHelper>(quoteHandles[Ndepos + i], (i + 1), (i + 7), 2,
                TARGET(), ModifiedFollowing, false, Actual360()));
    }

    // swaps
    auto euribor6m = ext::make_shared<Euribor>(6 * Months);
    for (Size i = 0; i < maximumMaturity; ++i) {
        instruments.push_back(ext::make_shared<SwapRateHelper>(
            quoteHandles[Ndepos + Nfra + i], (i + 1) * Years, TARGET(), Annual, ModifiedFollowing,
            Thirty360(Thirty360::European), euribor6m));
    }

    // build a piecewise yield curve
    using CurveType = PiecewiseYieldCurve<ZeroYield, Linear, IterativeBootstrap>;
    return Handle<YieldTermStructure>(
        ext::make_shared<CurveType>(referenceDate, instruments, Actual365Fixed()));
}

// creates the Swap portfolio, given the curve
std::vector<ext::shared_ptr<VanillaSwap>> setupPortfolio(Size portfolioSize, Size maximumMaturity, Handle<YieldTermStructure> curveHandle) 
{
    auto euribor6mYts = ext::make_shared<Euribor>(6 * Months, curveHandle);

    // set up a vanilla swap portfolio
    euribor6mYts->addFixing(Date(2, October, 2014), 0.0040);
    euribor6mYts->addFixing(Date(3, October, 2014), 0.0040);
    euribor6mYts->addFixing(Date(6, October, 2014), 0.0040);

    std::vector<ext::shared_ptr<VanillaSwap>> portfolio;
    MersenneTwisterUniformRng mt(42);

    for (Size j = 0; j < portfolioSize; ++j) 
    {
        Real fixedRate = mt.nextReal() * 0.10;
        Date effective(6, October, 2014);
        Date termination = TARGET().advance(effective, static_cast<Size>(mt.nextReal() * static_cast<double>(maximumMaturity) + 1.) * Years);
        Schedule fixedSchedule(effective, termination, 1 * Years, TARGET(), ModifiedFollowing, Following, DateGeneration::Backward, false);
        Schedule floatSchedule(effective, termination, 6 * Months, TARGET(), ModifiedFollowing, Following, DateGeneration::Backward, false);

        portfolio.push_back(ext::make_shared<VanillaSwap>(
            VanillaSwap::Receiver, 10000000.0 / portfolioSize, fixedSchedule, fixedRate,
            Thirty360(Thirty360::European), floatSchedule, euribor6mYts, 0.0, Actual360()));
    }

    return portfolio;
}


// create tape
using tape_type = Real::tape_type;
tape_type tape;

// price with sensitivities using AAD
double calculate_npv(Handle<YieldTermStructure> curveHandle, const std::vector<Real>& marketQuotes, std::vector<ext::shared_ptr<VanillaSwap>> portfolio, std::vector<double>& gradient, int use_aad)
{
    gradient.clear();

    // Loop over each swap and price it.
    auto pricingEngine = ext::make_shared<DiscountingSwapEngine>(curveHandle);
    double v = 0.0;
    for (auto& swap : portfolio) {
        swap->setPricingEngine(pricingEngine);

        auto npv = swap->NPV();
        if (use_aad == 1)
        {
            derivative(npv) = 1.0;
            tape.computeAdjoints();
            std::transform(marketQuotes.begin(), marketQuotes.end(), std::back_inserter(gradient), [](const Real& q) { return derivative(q); });
        }

		v += value(npv);
    }

    return v;
}

void SaveResults(const std::map<Size, std::tuple<double, double, double, double>>& results)
{
    std::filesystem::path dataPath("C:/Working/sharpe-src-v138/Tests/AADTests/data");
    std::ofstream sensitivities_result(dataPath / "sensitivities_result.csv");
   
	sensitivities_result << "PortfolioSize;npv_no_aad;time_no_aad;npv_aad;time_aad\n";
    for (const auto& [pSize, data] : results) 
    {
        sensitivities_result << pSize << ";" << std::get<0>(data) << ";" << std::get<1>(data) << ";" << std::get<2>(data) << ";" << std::get<3>(data) << "\n";
    }
    sensitivities_result.close();

}

int main() {
    try {
        std::vector<Size> portfolioSize;
        for(Size i=500; i<=10000; i+=500)
			portfolioSize.push_back(i);
        
        Size maxMaturity = 40;

        // market Quotes
        std::vector<double> marketQuotes;
        prepareQuotes(maxMaturity, marketQuotes);
        // convert double market quotes into AD type (Real is active double type)
        std::vector<Real> marketQuotesAD(marketQuotes.begin(), marketQuotes.end());

        // evaluation date
        Date referenceDate(2, January, 2015);
        Settings::instance().evaluationDate() = referenceDate;

        // Construct and price larger and larger portfolios
        // We make a map of NPV per portfoliosize
		std::map<Size, std::tuple<double, double, double, double>> results;

        for (int use_aad = 0; use_aad < 2; use_aad++)
        {
            for (auto pSize : portfolioSize)
            {
                if (use_aad == 1)
                {
                    // register independent variables with the tape as inputs and start a new recording
                    tape.clearAll();
                    tape.registerInputs(marketQuotesAD);
                    tape.newRecording();
                }
                // build curve. We rebuil in each loop.
                auto curveHandle = bootstrapCurve(Settings::instance().evaluationDate(), marketQuotesAD, maxMaturity);
               
                std::vector<double> gradient;
                auto start = std::chrono::high_resolution_clock::now();
                auto portfolio = setupPortfolio(pSize, maxMaturity, curveHandle);
                auto v = calculate_npv(curveHandle, marketQuotesAD, portfolio, gradient, use_aad);
                auto end = std::chrono::high_resolution_clock::now();
                auto time_sensi = static_cast<double>(std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()) * 1e-3;
             
                if(use_aad == 1 && pSize % 1000 == 0)
				    std::cout << "Time pricing " << pSize << " swaps: "<< time_sensi << "ms\n";
            
                auto it = results.find(pSize);
                if (it == results.end()) 
                {
					std::tuple<double, double, double, double> data = std::make_tuple(v, time_sensi, 0.0, 0.0);
                    results[pSize] = data;; 
                }
                else 
                {
					auto data = results[pSize];
					data = std::make_tuple(std::get<0>(data), std::get<1>(data), v, time_sensi);
					results[pSize] = data;
                }
            }
        }
       
        SaveResults(results);
        return 0;
    }
    catch (std::exception& e) {
        std::cerr << e.what() << std::endl;
        return 1;
    }
    catch (...) {
        std::cerr << "unknown error" << std::endl;
        return 1;
    }
}

[auto-differentiation-dev]

Thanks for the thorough write-up and for stress-testing AAD in QuantLib — very useful.

The scaling behaviour you're seeing isn't what we'd expect in a standard AAD workflow.

What's happening

The computeAdjoints call rolls back all operations on the tape up to the last newRecording. Because the tape is linear, the operations of each new swap pricing get appended to the tape, on top of the previous one, and computeAdjoints rolls back all the previous ones as well. Also, since derivatives aren't reset between swaps, they accumulate — leading to incorrect results (you'll actually get the cumulative sum of derivatives over all previous swaps).

Ways to address this

1. Use computeAdjointsTo
   Roll back only to a specific point, i.e. where the pricing of each swap starts.
   Track the position with getPosition and pass it to computeAdjointsTo.
   Also call clearDerivatives between swaps to avoid accumulation.
   Though due to QuantLib's lazy evaluation which may re-use cached results, care must be taken.

2. Reset the tape between swaps
   Use resetTo so only one swap NPV remains on the tape.
   Also call clearDerivatives between swaps to avoid accumulation.
   Note: QuantLib’s lazy evaluation may reuse cached results, so handle with care. If there are issues, perhaps reset to after the first was priced.

3. Fully reset the tape
   After each swap, reset the tape completely and start with a new newRecording.
   You'll need to register your inputs again and construct each swap individually.

We recommend option 2 as the most efficient approach.
Please let us know how it works for you.

[fredrik-fyring]

Thanks for the suggestions. I will try them out and report back.

[xcelerit-squad]

Hi @fredrik-fyring, do you have an update on this? Let us know if we can close this issue or if you need more help. Thanks

[auto-differentiation-dev]

Closing this because of inactivity.