finmath-lib, Strata, JQuantLib vs. NablaTensor

The last two articles measured one thing against itself. The compiled tape against the tape interpreter; the draw cache against no draw cache. Useful, and also the easiest kind of benchmark to win, because the author controls both sides.

So this one puts the same two tapes and the same seven-point spot ladder in front of the three Java libraries a quant would reach for instead of NablaTensor, and runs each of them the way its own documentation says to.

• finmath-lib 6.1.9 — Christian Fries' library, the standard Java mathematical-finance stack for Monte Carlo, LIBOR market models and xVA. It has backward-mode automatic differentiation built in — a separate extension module until it was folded into the core library at version 3.3 — so it gets to use it.
• Strata 2.12.74 — OpenGamma's Apache-2.0 market-risk library, used in production at banks for curves, pricing and sensitivities.
• JQuantLib — the Java port of QuantLib, in the 0.3.0 repackaging that is on Maven Central.

The rules, stated before the numbers so that they cannot be adjusted after:

1. Every library runs as it ships. If it offers adjoint Greeks, it uses them. If it memoises random draws, that is used too. If it exposes a thread count, it gets eight threads like the NablaTensor engines. If it does not, it runs on one, and the table says so — I did not wrap anything in an executor to manufacture an eight-thread number. What a library can do is the measurement; what I could build on top of it is not.
2. Same product, same market, same schedule: S = K = 100, σ = 20%, r = 3%, T = 1, the one-step European call and the 252-fixing arithmetic Asian call, log-Euler on a uniform grid.
3. Same ladder: spots 98, 99, 99.5, 100, 100.5, 101, 102, one seed, warm-up on a different seed, three ladders, median of the six warm calls. The NablaTensor rows are the ones already published in the draw-cache article, plus a one-thread cpu-jit row measured for this piece so that there is a per-core comparison in both directions.

Everything below was measured on referential machine — Ryzen 7 8845HS, Zulu OpenJDK 25.0.1 — with the three libraries pulled from Maven Central at the versions above.

What each library turned out to be able to do

This section was supposed to be one paragraph. It is the most important section in the article.

finmath-lib can do all of it. BlackScholesModel plus EulerSchemeFromProcessModel plus a payoff written against its RandomVariable type is the Asian in nine lines. Wrap the model's inputs in RandomVariableDifferentiableAADFactory and the same nine lines produce delta, strike sensitivity, vega and rho from one getGradient() call — a real backward sweep over an operator graph, the same idea as NablaTensor's tape. And BrownianMotionFromMersenneRandomNumbers generates its increments lazily and keeps them, so two models built on the same BrownianMotion object share draws: finmath has had a draw cache since before NablaTensor existed, and it is on by default.

What finmath does not have, on this workload, is a second core. EulerSchemeFromProcessModel does have multithreading — a system property, isUseMultiThreadding, on by default — but it parallelises across model components, and a Black-Scholes model has one. The arithmetic on a RandomVariableFromDoubleArray is a plain sequential loop. Measured with the serial collector so that GC threads could not confuse the reading: 8.25 s of user time in 9.91 s of wall clock. One thread, as shipped. (An eight-way path shard is about forty lines of your own code and works fine; I wrote it, it is not in the tables, because it is my code and not the library's.)

Strata has no Monte-Carlo engine. None. It is a market-risk library built on closed forms, curves and calibration, and its equity-option content is BlackScholesFormulaRepository: price, delta, dual delta, gamma, vega, rho, theta and a dozen higher-order Greeks, each a static method on six doubles and a call/put flag. For the European that is the right tool and it appears below. For the Asian there is nothing to call.

JQuantLib has an MCEuropeanEngine class, and I was looking forward to it. Its constructor, in the 0.3.0 jar:

public MCEuropeanEngine();
  Code:
     4: new           #2   // class java/lang/UnsupportedOperationException
     8: ldc           #3   // String work in progress
    13: athrow

MCVanillaEngine has an empty constructor and nothing else; the Asian engines present are the two closed-form geometric-average ones. The Monte-Carlo half of the QuantLib port was never finished. AnalyticEuropeanEngine does work, and returns QuantLib's full Greek set from one BlackCalculator pass, so JQuantLib also appears in the European table — as a closed form.

So the honest shape of the comparison is: on the European, five contenders, two of them not simulating anything; on the Asian, the one that actually needs Monte Carlo, one library besides NablaTensor can run it at all.

The same thing as a figure, with the two questions a reader asks before any benchmark — can I <dependency> it, and what hardware does it use — added:

The one-step European

Twenty million paths, fp64. finmath is measured once per configuration on one thread; the NablaTensor rows are from the draw-cache article, plus the one-thread cpu-jit rows measured for this one.

library	threads	Greeks	draws	price only	value + Greeks	peak RSS, price / Greeks
finmath-lib 6.1.9	1	AAD, 4	fresh BrownianMotion per call	33.7 Mpath/s	17.2	11.1 / 16.8 GB
finmath-lib 6.1.9	1	AAD, 4	one BrownianMotion, memoised	68.3	23.7	10.1 / 16.5 GB
cpu-jit	1	adjoint, 5	regenerated	23.8	15.1	67 / 70 MB
cpu-jit + DRAW_CACHE	1	adjoint, 5	cached	40.5	23.6	351 / 281 MB
cpu-jit	8	adjoint, 5	regenerated	155.6	98.9	68 / 72 MB
cpu-jit + DRAW_CACHE	8	adjoint, 5	cached	285.4	150.1	369 / 375 MB
simd	8	adjoint, 5	regenerated	340.9	250.7	127 / 553 MB
simd + DRAW_CACHE	8	adjoint, 5	cached	632.1	340.9	428 MB / 1.2 GB

The last column is the peak resident set of the whole process — every thread, the JVM, and the draw cache — for the price-only run and for the Greeks run, read with GNU time from a separate one-ladder run of each configuration. The NablaTensor rows ran with the default heap; the finmath rows with the -Xmx16g the Asian below needs, and finmath fills whatever heap it is given: every RandomVariable operation allocates a fresh 160 MB double[20,000,000], and G1 lets them pile up until the ceiling. The same European ladder completes in a 1 GB heap on the price pass and a 4 GB heap on the Greeks pass, so read the finmath cells as "what it takes when it can", not as a floor. (A spot check worth its own article: with -Xmx1g the price-only pass ran 2.8× faster than with -Xmx16g on this box, because the collector reuses the same few regions instead of touching fresh pages for every array. The 68.3 above was measured with a 14 GB heap.)

And the two closed forms, which do not simulate and are listed in their own unit because a "path" is not a thing they have:

library	threads	what a valuation is	valuations / s	peak RSS
Strata 2.12.74 BlackScholesFormulaRepository	1	price + delta, dual delta, vega, rho, theta (six static calls)	4.10 M	60 MB
JQuantLib 0.3.0 AnalyticEuropeanEngine	1	NPV() + delta, vega, rho, theta from one BlackCalculator	4.08 M	527 MB

Three things to take from this, and the first is the one I did not expect.

On one core, finmath is as fast as the compiled tape, and on the price pass it is faster. 33.7 against 23.8 million paths per second with fresh draws; 68.3 against 40.5 with cached ones. On the Greeks pass, memoised draws on both sides, it is a dead heat: 23.7 against 23.6. NablaTensor's whole first article was about the cost of interpreting a tape per path, and here is a library that does not compile anything, does not generate bytecode, and matches it.

The reason is that finmath does not interpret per path either. Its RandomVariable is an array of paths, and every operation — mult, exp, sub, floor — is one tight loop over twenty million doubles. That is the NumPy design: dispatch once per operation, not once per operation per path, so the dispatch cost that cpu-jit had to compile away is amortised over the whole array instead. On a one-step tape with a dozen operations, each of them a streaming pass over 160 MB, it is a perfectly good way to run a Monte Carlo. The draw generation is also cheaper: Mersenne Twister plus an inverse normal CDF, against Philox's ten rounds and Box-Muller's log, sqrt and sincos.

The memoised draws are worth the same thing in both libraries. finmath's price pass goes 2.0× faster on cached increments, its Greeks pass 1.4×; the draw-cache article measured 1.8× and 1.5× for cpu-jit. That is not a coincidence, it is Amdahl on the same RNG share, and it is a quiet confirmation that the two libraries are doing comparable amounts of work per path once the random numbers are out of the way.

The 6.3× between them at the end of the table is the thread count and nothing else. cpu-jit on eight threads with the cache does 150.1 on the Greeks pass; finmath, at parity per core, does 23.7. Every bit of that gap is cores finmath cannot use, and every bit of simd's further 2.3× is the Vector API. On this tape the compiled kernel wins by being parallel, not by being compiled.

And the closed forms: 4.1 million complete valuations per second, price and five Greeks, on one thread. That is 240 nanoseconds each. Twenty million Monte-Carlo paths at 632 million per second take 32 milliseconds to reach a price with a standard error of 0.003; Strata reaches the exact answer in a quarter of a microsecond. Nobody prices a vanilla European by simulation, and the row is here to say so in numbers. The Monte-Carlo rows exist because of the next table.

The 252-step Asian

The tape that needs simulation: 1,536 nodes, 252 draws per path, an arithmetic average that no closed form covers. Strata has nothing to call. JQuantLib's engine throws. finmath runs it — at 200,000 paths rather than the 300,000 of the NablaTensor rows, because its AAD operator graph keeps a double[paths] per node and 300,000 paths did not fit in a 16 GB heap. Throughput is per path and does not depend on the count at this size, but the difference is stated.

library	threads	Greeks	draws	price only	value + Greeks	peak RSS, price / Greeks
finmath-lib 6.1.9	1	AAD, 4	fresh BrownianMotion per call	0.125 Mpath/s	0.089	4.7 / 7.2 GB
finmath-lib 6.1.9	1	AAD, 4	one BrownianMotion, memoised	0.338	0.174	4.8 / 7.7 GB
cpu-jit	1	adjoint, 5	regenerated	0.26	0.22	76 / 112 MB
cpu-jit + DRAW_CACHE	1	adjoint, 5	cached	0.84	0.51	995 / 829 MB
cpu-jit	8	adjoint, 5	regenerated	1.69	1.44	80 / 119 MB
cpu-jit + DRAW_CACHE	8	adjoint, 5	cached	4.62	3.07	1.2 / 1.3 GB
simd	8	adjoint, 5	regenerated	4.29	3.12	191 / 419 MB
simd + DRAW_CACHE	8	adjoint, 5	cached	7.69	3.97	1.4 / 1.4 GB
Strata 2.12.74	—	—	—	no Monte-Carlo engine
JQuantLib 0.3.0	—	—	—	MCEuropeanEngine: "work in progress"

Same memory column as above, whole process, price run / Greeks run. The NablaTensor rows are at 300,000 paths, so their draw cache is 300,000 × 252 doubles, 605 MB, and that is what separates the cached rows from the uncached ones. The finmath Greeks pass at 200,000 paths is 7.7 GB and does not fit at 300,000 paths with a 16 GB heap; that is not the heap filling up for lack of collection, it is the operator graph, as the next section says.

Read the two bold columns. On the Greeks pass, best configuration against best configuration, cpu-jit on eight threads is 17.6× finmath. But the one-thread cpu-jit row is there so that the ratio can be taken apart: cpu-jit on one core with the cache does 0.51 against finmath's 0.174 — a 2.9× per-core gap — and eight threads then multiply it by 6.0. On the European the per-core gap was 1.0. Something about the long tape costs finmath a factor of three per core that the short tape did not, and the memory column says what.

Where the three-per-core goes

The array-of-paths design has a cost that the one-step tape hid and the 252-step tape exposes: every intermediate is materialised. A RandomVariable operation reads one or two arrays of 200,000 doubles from memory and writes a new one. The Asian is 252 Euler steps of drift, diffusion, exp, multiply and running sum, then the average and the payoff — on the order of 1,500 vector operations, each a 1.6 MB write. The EulerSchemeFromProcessModel keeps every step's state (it has to, the payoff reads all 252 fixings), and the AAD factory keeps every operand of every operation, because the backward sweep will need them. That is the 7.7 GB in the last column: not a leak, the operator graph doing exactly what an operator graph over arrays has to do.

The compiled tape holds the same 1,536 values per path — 12 KB, in L1 — and throws them away when the path is done. It touches main memory for the draws and for nothing else. That is the difference between 119 MB and 7.7 GB, and it is also, once the working set is bigger than any cache on the chip, the difference between 0.51 and 0.174 million paths per second on a single core. Same Wengert list, same reverse sweep, same Baur–Strassen bound on the adjoint cost — finmath's Greeks pass is 1.9× its price pass with memoised draws, cpu-jit's is 1.66× — but one of them is transposed, and the transposed one streams its tape through DRAM.

The draws, again

finmath's memoised BrownianMotion lifts its Asian price pass 2.7× and its Greeks pass 2.0×. The draw-cache article measured 2.7× and 2.1× for cpu-jit. Two libraries with different generators — Mersenne Twister with an inverse CDF against Philox with Box-Muller — different array layouts and no shared code arrive at the same lever with the same leverage, because on a 252-draw tape the random numbers are around two thirds of a price path whichever way you compute the rest. If there is one number in this article that generalises past the two libraries measured, it is that one.

Did they agree?

Different generators mean different noise, so nothing here is bit-comparable across libraries; what should agree is the estimate within its own standard error, and the Greeks within theirs.

	price	delta	dV/dK	vega	rho	dV/dT
European, closed form (Strata, JQuantLib)	9.4134	0.59871	−0.50457	38.667	50.457	5.380
European, finmath, 20 M paths, AAD	9.4188	0.59883	−0.50464	38.697	50.464	—
European, cpu-jit, 20 M paths, adjoint	9.4122	0.59870	−0.50458	38.661	50.458	5.380
Asian, finmath, 200 k paths, AAD	5.3140	0.56210	−0.50896	22.465	23.603	—
Asian, cpu-jit, 300 k paths, adjoint	5.3010	0.56163	−0.50862	22.386	23.591	2.946

(dV/dT is the derivative with respect to maturity; Strata and QuantLib report its negative and call it theta, −5.380.)

The European standard error at twenty million paths is 0.0032 as reported by the engine; finmath sits 1.7 of them above the closed form and cpu-jit 0.4 below, which is what two independent estimates of the same number look like. The Asian's standard error is 0.014 at 300,000 paths, and the two estimates are 0.013 apart. The Greeks agree to three figures everywhere they overlap; the one that does not overlap is dV/dT, because finmath's time grid is not something its AAD can differentiate through, while for the tape maturity is just another input node.

What to actually use

The ranking of the things that decided this comparison, in the order they decided it:

1. Whether the library can run the product at all. Two of three could not run the Asian. No amount of per-core efficiency matters to a row that says "work in progress".
2. Threads. 6× of the 17.6× on the Asian, all of the 6.3× on the European. finmath is a single-threaded library on a single-component model, as shipped; you can shard it yourself, and then the number you get is a measure of your sharding.
3. Per-path versus per-array. A wash on a twenty-six-node tape; 2.9× per core and 7.7 GB against 119 MB on a 1,536-node one. This is the only line in the list that is about NablaTensor's design rather than its configuration, and it only shows up on the tape that needed simulation in the first place.
4. The draw cache. 2–2.7× on the price pass in both libraries, and finmath's is on by default. If you take one engineering idea from this series to a codebase that is not NablaTensor, take that one.
5. Adjoint Greeks. Not a differentiator between finmath and NablaTensor, because both have them and both pay about 1.7–1.9× a price pass for every Greek at once. A differentiator against everything that bumps.

If your book is closed-form products and curves, Strata is the production library and this article has nothing to tell you; 4.1 million valuations a second with the Greeks attached is the correct way to price a European. If your book is path-dependent and single-threaded is acceptable — or you are content to shard it yourself — finmath is a serious, honest, well-designed library that has had adjoint Greeks and cached draws for years, and it will tie the compiled tape per core on anything short. If your book is path-dependent, your tapes are long and you have cores, the tape kernel is 18× on this box before anyone touches a GPU.

Try it

The harnesses are plain main methods against Maven Central artefacts — no build of any of the four libraries is required:

<dependency><groupId>net.finmath</groupId><artifactId>finmath-lib</artifactId><version>6.1.9</version></dependency>
<dependency><groupId>com.opengamma.strata</groupId><artifactId>strata-pricer</artifactId><version>2.12.74</version></dependency>
<dependency><groupId>co.dv01.jquantlib</groupId><artifactId>core</artifactId><version>0.3.0</version></dependency>
<dependency><groupId>com.nablatensor</groupId><artifactId>nablatensor-quant</artifactId><version>0.1.0</version></dependency>
<dependency><groupId>com.nablatensor</groupId><artifactId>nablatensor-simd</artifactId><version>0.1.0</version></dependency>

nablatensor-quant pulls in cpu-jit; simd is the second artefact and needs --add-modules jdk.incubator.vector.

finmath's Asian with adjoint Greeks, in full — this is the code that produced its rows, minus the timing loop:

var f   = new RandomVariableDifferentiableAADFactory(new RandomVariableFromArrayFactory(true));
var rvS = f.createRandomVariable(spot);  var rvK = f.createRandomVariable(100.0);
var rvV = f.createRandomVariable(0.20);  var rvR = f.createRandomVariable(0.03);
var td  = new TimeDiscretizationFromArray(0.0, 252, 1.0 / 252);
var bm  = new BrownianMotionFromMersenneRandomNumbers(td, 1, paths, seed, new RandomVariableFromArrayFactory(true));
var model = new BlackScholesModel(rvS, rvR, rvV, f);
var sim   = new MonteCarloAssetModel(model, new EulerSchemeFromProcessModel(model, bm));

RandomVariable sum = sim.getAssetValue(td.getTime(1), 0);
for (int i = 2; i <= 252; i++) sum = sum.add(sim.getAssetValue(td.getTime(i), 0));
RandomVariable value = sum.div(252).sub(rvK).floor(0.0)
    .div(sim.getNumeraire(1.0)).mult(sim.getNumeraire(0.0));

Map<Long, RandomVariable> g = ((RandomVariableDifferentiable) value).getGradient();
double price = value.getAverage(), delta = g.get(rvS.getID()).getAverage();   // etc.

Reuse bm across the seven spots and you have finmath's draw cache; build a new one per spot and you have the "fresh" row. Set -Xmx16g for 200,000 paths and watch the resident set with /usr/bin/time -v, because the number in the last column of the Asian table is the one that will decide whether your ladder fits on your box before the throughput does.

Source: bench/three-java-quant-libraries-one-ladder (the three harnesses and their pom.xml, on the main engine repo), JitReplay.java (the NablaTensor side), and for finmath RandomVariableDifferentiableAADFactory, BrownianMotionFromMersenneRandomNumbers and EulerSchemeFromProcessModel (the memoised increments and the per-component executor).