life_modeling_problem.md

library(data.table)

q <- c(0.001,
       0.002,
       0.003,
       0.003,
       0.004,
       0.004,
       0.005,
       0.007,
       0.009,
       0.011)

w <- c(0.05,
       0.07,
       0.08,
       0.10,
       0.14,
       0.20,
       0.20,
       0.20,
       0.10,
       0.04)

P <- 100
S <- 25000
r <- 0.02

dt <- as.data.table(cbind(q,w))

npv <- function(cf, r, S, P) {
  cf[, inforce := shift(cumprod(1 - q - w), fill = 1)
  ][, lapses := inforce * w
  ][, deaths := inforce * q
  ][, claims := deaths * S
  ][, premiums := inforce * P
  ][, ncf := premiums - claims
  ][, d := (1/(1+r))^(.I)
  ][, sum(ncf*d)]
}

npv(dt,r,S,P)
#> [1] 50.32483

microbenchmark::microbenchmark(npv(dt,r,S,P))
#> Unit: milliseconds
#>              expr    min      lq     mean  median      uq    max neval
#>  npv(dt, r, S, P) 2.5791 2.71035 2.964293 2.85625 3.10385 6.0357   100

^{Created on 2021-01-15 by the reprex package (v0.3.0)}

@alecloudenback, thoughts on reducing allocations? StaticArrays?

StaticArrays help when the size of the array is less than about 15 items, IIRC. Since the sample here is unrealistically short (e.g. q is often much longer). So it would probably help here, but would be misleading for more "real world" problem sizes.

Here's a small tweak to the code, which avoids the more expensive concatenation ([arr1;arr2]) and modifies the same array that is generated with the cumprod:

function npv2(q,w,P,S,r) 
	inforce = append!([1.0],cumprod(1 .- q .- w)[1:end-1])
  	ncf = inforce .* P .- inforce .* q .* S
 	d = (1 ./ (1 + r)) .^ (1:length(ncf))
  	return sum(ncf .* d)
end

which is close to ~8x speedup:

BenchmarkTools.Trial: 
  memory estimate:  1.12 KiB
  allocs estimate:  8
  --------------
  minimum time:     413.945 ns (0.00% GC)
  median time:      427.553 ns (0.00% GC)
  mean time:        485.135 ns (9.34% GC)
  maximum time:     15.881 μs (95.88% GC)
  --------------
  samples:          10000
  evals/sample:     199

And inserting an @view where the cumprod array is sliced makes the operation non-allocating and eeks out a touch more performance:

function npv2(q,w,P,S,r) 
	inforce = append!([1.0],@view cumprod(1 .- q .- w)[1:end-1])
  	ncf = inforce .* P .- inforce .* q .* S
 	d = (1 ./ (1 + r)) .^ (1:length(ncf))
  	return sum(ncf .* d)
end

BenchmarkTools.Trial: 
  memory estimate:  992 bytes
  allocs estimate:  7
  --------------
  minimum time:     390.470 ns (0.00% GC)
  median time:      398.723 ns (0.00% GC)
  mean time:        448.883 ns (8.84% GC)
  maximum time:     13.341 μs (95.96% GC)
  --------------
  samples:          10000
  evals/sample:     202

Also - just realised we can't compare results as we're not all using the same hardware!

@lewisfogden - I'd say that's partly true. In the cases here, we are comparing orders of magnitude difference. The relative performance of modern processors is noise on this scale:

code	relative	absolute timing (mean)
Julia	1x	500 nanoseconds
Python numpy	42x	21 microseconds
R vectorized	140x	70 microseconds
R data.table	6000x	3 milliseconds

Notes: using mean times, since that's more relevant for overall time on larger problems (you'd want to know how long it would take to complete the combined problem, not just the typical (median) sub-problem). Also not comparing to the memoized version posted above since that's a different technique.

I'm on an M1 Macbook Air (where Julia has to run on Rosetta2 because there's not a native version available. For calibration's sake, I ran the "R vectorized" version posted above on my machine (running on Rosetta2 as well):

Unit: microseconds
                      expr   min   lq     mean median     uq      max neval
 base_r_npv(q, w, P, S, r) 7.042 7.25 68.62107  7.418 9.6045 5989.668   100

The median is pretty close to @Houstonwp above. I don't understand why the mean is so much higher in our benchmarks than the median? I still think if one is trying to figure out the best solution to crunching a problem at scale, the mean is a better representation... but what's going on here?

I didn't try running the python version since the python binary and environment management is so notoriously tricky... I'm waiting for an ARM-native homebrew, which seems to be the best practice of managing python on Mac

I don't understand why the mean is so much higher in our benchmarks than the median

I've been running them in a Ubuntu VM on Windows 10 (using WSL) - the 10 manual runs I did for nim indicated that there was a spike every so often (e.g. some background operation running). The mean is skewed by that high 5989 maximum - it is contributing ~5989/100 = 60 to the mean, which would be around 8 microsecs otherwise.

If you are doing homogeneous operations then the median is probably a more useful indicator purely for benchmarking - if it was at scale and you wanted to test a variety of policies then agree mean (or a truncated mean) would be better to indicate your production run time.

I had a look at the Julia Benchmarks - they roughly tie up to your table. LuaJIT (which I've never heard of) gives Julia a good run for its money on the JIT front. Downside on these is that its probably harder to interface them into Python if you are building a production pipeline.

Here is the Rust implementation for reference:

pub fn npv(mortality_rates: &Vec<f64>, lapse_rates: &Vec<f64>, interest_rate: f64, sum_assured: f64, premium: f64, init_pols: f64, term: Option<usize>) -> f64 {

    let term = term.unwrap_or_else(|| mortality_rates.len());
    let mut result = 0.0;
    let mut inforce = init_pols;
    let v: f64 = 1.0 / (1.0 + interest_rate);

    for (t, (q, w)) in mortality_rates.iter().zip(lapse_rates).enumerate() {
        let no_deaths = if t < term {inforce * q} else {0.0};
        let no_lapses = if t < term {inforce * w} else {0.0};
        let premiums = inforce * premium;
        let claims = no_deaths * sum_assured;
        let net_cashflow = premiums - claims;
        result += net_cashflow * v.powi(t as i32);
        inforce = inforce - no_deaths - no_lapses;
    }

    result
}

And the benchmark:

npv                     time:   [362.23 ns 364.30 ns 366.35 ns]
Found 9 outliers among 100 measurements (9.00%)
  7 (7.00%) low mild
  1 (1.00%) high mild
  1 (1.00%) high severe

The reason there are three measurements above is that it's showing the upper and lower ranges. The Rust benchmarking system is a statistical one.

I wrote this fast so I may have made a mistake. I also did not try to optimize it so there's some performance gain to be had for sure. Even though it looks serial I'm sure the compiler is auto-vectorizing it.

Nice, @paddyhoran. I have not used rust so played with it here: https://play.rust-lang.org/?version=stable&mode=debug&edition=2018&gist=255a33e84c70e3bcd43258ddbaa66f27 Looks like the answer matches the other above.

Also looks like your algorithm for the problem is more efficient overall. It's over 2x faster than the prior Julia fastest Julia version above. Code and benchmark:

function npv3(q,w,P,S,r,term=nothing)
	term = isnothing(term) ? length(q) + 1 : term + 1
	inforce = 1.0
	result = 0.0
	
	for (t,(q,w)) in enumerate(zip(q,w))
		deaths = t < term ? inforce * q : 0.0
		lapses = t < term ? inforce * w : 0.0
		premiums = inforce * P
		claims = deaths * S
		ncf = premiums - claims
		result += ncf * (1 / ( 1 + r)) ^ (t-1)
		inforce = inforce - deaths - lapses
	end
	
  	return result
end

BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     135.023 ns (0.00% GC)
  median time:      135.262 ns (0.00% GC)
  mean time:        137.513 ns (0.00% GC)
  maximum time:     1.327 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     873

Updated Rankings for the "Life Modeling Problem":

code	algorithm	relative	absolute timing (mean)
Julia	accumulating for loop	1x	135 nanoseconds
Rust	accumulating for loop	3x	360 nanoseconds
Julia	vector	4x	500 nanoseconds
Python numpy	vector	155x	21 microseconds
R vectorized	vector	520x	70 microseconds
R data.table	data.table vector	22000x	3 milliseconds