14Accelerators for Quantitative Work

A bank moved its overnight counterparty-risk Monte Carlo to graphics processors and it ran many times faster. The same bank tried the same devices for the intraday pricing of single trades, and nothing moved: each request spent longer crossing to the device and back than computing on it. Both results were predictable from three numbers per workload — how much arithmetic it does, how many bytes it moves, and how many bytes must cross to the device — and two numbers per machine. This chapter measures those numbers on the laptop’s processor, takes a data-centre device’s from its product page, and predicts both outcomes before anyone buys anything.

14.1 What a graphics processor is good at

A graphics processor (an accelerator in the sense of Book 12, chapter 23) has thousands of simple arithmetic units running the same instruction on different data, and memory attached to it that is several times faster than a server’s. It is good at the same thing SIMD instructions are good at on a CPU (Book 13, chapter 14), at a much larger scale: a long loop of identical, independent arithmetic on arrays already in its memory. It is bad at three things quantitative code often does.

The third is size. A computation that takes a microsecond on a CPU cannot be faster on a device whose launch and synchronisation alone take several: small work, like a single trade’s price, pays a fixed cost that large work amortises.

14.2 The roofline model

The roofline turns a question about hardware into two numbers about code. The chapter’s model (firm.roofline) adds the transfer and the fixed cost: on a device, a kernel takes at least $L + X/\ell + \max(F/P, M/B)$ seconds, with $F$ its operations, $M$ its memory bytes, $X$ the bytes that cross a link of bandwidth $\ell$, and $L$ the launch and synchronisation latency of a request.

@dataclass(frozen=True)
class Roof:
    name: str
    peak: float                       # flop/s
    bandwidth: float                  # byte/s
    link: float = math.inf            # host link, byte/s
    launch: float = 0.0               # s per request (launch and synchronisation)
    source: str = ""

    @property
    def ridge(self) -> float:
        return self.peak / self.bandwidth

    def attainable(self, intensity: float) -> float:
        return min(self.peak, intensity * self.bandwidth)

    def time(self, k: Kernel) -> float:
        return max(k.flops / self.peak, k.bytes / self.bandwidth)

    def end_to_end(self, k: Kernel) -> float:
        return self.launch + k.transfer / self.link + self.time(k)

14.3 Measuring the laptop’s roof

One core of this laptop (an Intel Core Ultra 7 155H, under WSL2) was measured: the peak by a $2000\times2000$ double-precision matrix product on one BLAS thread, the bandwidth by a STREAM-style triad $a_i = b_i + s c_i$ in C++20 on arrays of 128 MB, counting 24 bytes per element. The core reaches 66.9 GFLOP/s and 24.0 GB/s: its ridge is at 2.79 operations per byte.

    std::vector<double> a(n, 0.0), b(n, 1.0), c(n, 2.0);
    const double s = 3.0;
    double best = 1e30;
    for (int r = 0; r < reps; ++r) {
        auto t0 = std::chrono::steady_clock::now();
        for (std::size_t i = 0; i < n; ++i) a[i] = b[i] + s * c[i];
        auto t1 = std::chrono::steady_clock::now();
        double dt = std::chrono::duration<double>(t1 - t0).count();
        if (dt < best) best = dt;
        b[r % n] = a[(r * 7) % n];           // keep the loop from being optimised away
    }
    std::printf("%.3f %.9f\n", 24.0 * static_cast<double>(n) / best / 1e9, a[n / 2]);

Four reference kernels were timed on the same core (Table 14.1): a step of geometric Brownian motion over $2^{23}$ paths ($S \mathrel{*}= e^{a + bZ}$), a basket payoff over 500 000 paths of 20 assets, the aggregation of a cube of trade values (2 000 paths, 50 dates, 200 trades) into an expected-positive-exposure profile, and the matrix product. The operations are counted by the convention of firm.roofline: additions and multiplications count one, and so does a call to exp, although it costs many cycles.

Three of the four are memory-bound, and the two that are simple loops reach 81 and 87% of the bandwidth roof: they are done, and only more bandwidth would speed them up. The path step sits far below its roof because the roofline counts operations and exp is one operation that takes many cycles; its true limit is the rate of the exponential. The matrix product, the only compute-bound kernel, defines the peak: it reaches 96% of it here only because the peak is its best over more runs than the table’s five.

14.4 Monte Carlo and risk on a device

The overnight workload is an exposure Monte Carlo for one netting set: 10 000 paths of 50 risk factors on 120 monthly dates, 5 000 trades valued on every path and date (ten operations each: a linear trade on one factor, discounted), and the values aggregated into an expected-exposure profile (Book 6, chapters 17 and 18). It performs 66 billion operations and moves 98 GB through memory: memory-bound on both machines. Only the trades’ descriptions go to the device (100 bytes each) and only the profile comes back; the paths are generated where they are used, with a counter-based generator (Book 4, chapter 26) so that the device’s paths are the CPU’s.

def exposure_job() -> list[R.Kernel]:
    n = PATHS * DATES * TRADES
    factors = PATHS * DATES * FACTORS
    paths = R.k_path(factors)
    value = R.Kernel("valuation", VALUE_FLOPS * n, 8.0 * (n + factors), TRADES * TRADE_BYTES)
    agg = R.k_exposure(PATHS, DATES, TRADES)
    return [paths, value, agg]


def overnight(cpu: R.Roof, dev: R.Roof = DEV) -> dict:
    ks = exposure_job()
    t_cpu = sum(cpu.time(k) for k in ks)
    moved = sum(k.transfer for k in ks) + DATES * 8.0            # trades in, the profile out
    t_dev = dev.launch * len(ks) + moved / dev.link + sum(dev.time(k) for k in ks)
    s = t_cpu / t_dev
    return {"cpu_s": t_cpu, "dev_s": t_dev, "speedup": s, "amdahl": R.amdahl(DEVICE_SHARE, s),
            "bytes": sum(k.bytes for k in ks), "flops": sum(k.flops for k in ks)}


def request(cpu: R.Roof, n: int = 1, dev: R.Roof = DEV) -> dict:
    k = R.Kernel("price request", REQUEST_FLOPS * n, 0.8 * REQUEST_FLOPS * n, REQUEST_BYTES * n)
    t_cpu, t_dev = cpu.time(k), dev.end_to_end(k)
    return {"n": n, "cpu_s": t_cpu, "dev_s": t_dev, "speedup": t_cpu / t_dev}

On the core, at its roof, the job takes 4.09 seconds; on the device, 29.3 milliseconds including three launches and the transfers: 140 times faster, the ratio of the bandwidths, because the job is memory-bound. If 95% of the overnight process is this kind of work and the rest (reading positions, writing reports) stays on the host, Amdahl’s law caps the whole process’s gain at $1/(0.05 + 0.95/140) = 17.6$ times. Against a server socket of many cores the per-core factor shrinks: memory-bound work gains the ratio of the device’s bandwidth to the socket’s, not to one core’s.

14.5 Training, and what the machine-learning book already said

Book 12 (chapter 23) covered the accelerator’s home ground: training neural networks is dominated by matrix products, compute-bound with high intensity, run in mixed precision on tensor units, and split across devices by data parallelism. The same roofline explains it: a matrix product of size $n$ has intensity $n/12$, far above any ridge, and a device’s reduced-precision peak is many times its double-precision one. Quantitative finance’s Monte Carlo is usually the other kind of kernel — low intensity, double precision — and gains the bandwidth ratio, not the peak ratio.

14.6 Deciding with a cost model

A single intraday price request is small: ten thousand operations (a hundred cash flows on an interpolated curve) and a kilobyte crossing the link (the trade and its result). On the core it takes 0.33 microseconds; on the device, 10.0 — almost all of it the launch and synchronisation latency, which the chapter takes as 10 microseconds (an assumption, not a datasheet figure). The device is 30 times slower. Batched, the fixed cost is shared (Figure 14.3): the device breaks even at 31 trades per request and approaches, for very large batches, 32.7 times the core — a ceiling set by the link, because every trade’s kilobyte still crosses it.

14.7 Tutorial: two rooflines and a decision

14.8 Build: the roofline model

14.9 Exercises

14.10 Problem: Faster Overnight, Not Intraday

14.11 Interview questions