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.

PortfolioSize	time without XAD [ms]	time 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. Is this to be expected or could it be improved on somehow? I built the example below with "windows-xad clang-debug" and "windows-xad-clang-release" configuration on Windows 11.

Kind Regards
Fredrik

#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]

Closing this as it has been moved to in issue in #31