XIRR performance implrovements (faster deriv calculation)

CHANGELOG.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## not released
+
+- XIRR Performance improvements
+
 ## [0.10.2] - 2024-01-27
 
 - (X)IRR Performance improvements

README.md

@@ -32,7 +32,7 @@ and the [implementation from the Stack Overflow](https://stackoverflow.com/a/115
 
 ![bench](https://raw.githubusercontent.com/Anexen/pyxirr/main/docs/static/bench.png)
 
-PyXIRR is much faster in XIRR calculation than the other implementations.
+PyXIRR is much faster than the other implementations.
 
 Powered by [github-action-benchmark](https://github.com/benchmark-action/github-action-benchmark) and [plotly.js](https://github.com/plotly/plotly.js).
 
@@ -61,16 +61,16 @@ xirr(['2020-01-01', '2021-01-01'], [-1000, 1200])
 
 # Multiple IRR problem
 
-The Multiple IRR problem occur when the signs of cash flows change more than
+The Multiple IRR problem occurs when the signs of cash flows change more than
 once. In this case, we say that the project has non-conventional cash flows.
 This leads to situation, where it can have more the one IRR or have no IRR at all.
 
-PyXIRR's approach to the Multiple IRR problem:
+PyXIRR addresses the Multiple IRR problem as follows:
 
 1. It looks for positive result around 0.1 (the same as Excel with the default guess=0.1).
 2. If it can't find a result, it uses several other attempts and selects the lowest IRR to be conservative.
 
-Here is an example of how to find multiple IRRs:
+Here is an example illustrating how to identify multiple IRRs:
 
 ```python
 import numpy as np

src/core/optimize.rs

@@ -11,13 +11,42 @@ where
     let mut x = start;
 
     for _ in 0..MAX_ITERATIONS {
-        let res = f(x);
+        let y = f(x);
 
-        if res.abs() < MAX_ERROR {
+        if y.abs() < MAX_ERROR {
             return x;
         }
 
-        let delta = res / d(x);
+        let delta = y / d(x);
+
+        if delta.abs() < MAX_ERROR {
+            return x - delta;
+        }
+
+        x -= delta;
+    }
+
+    f64::NAN
+}
+
+// a slightly modified version that accepts a callback function that
+// calculates the result and the derivative at once
+pub fn newton_raphson_2<Func>(start: f64, fd: &Func) -> f64
+where
+    Func: Fn(f64) -> (f64, f64),
+{
+    // x[n + 1] = x[n] - f(x[n])/f'(x[n])
+
+    let mut x = start;
+
+    for _ in 0..MAX_ITERATIONS {
+        let (y0, y1) = fd(x);
+
+        if y0.abs() < MAX_ERROR {
+            return x;
+        }
+
+        let delta = y0 / y1;
 
         if delta.abs() < MAX_ERROR {
             return x - delta;

src/core/scheduled/xirr.rs

@@ -1,7 +1,7 @@
 use super::{year_fraction, DayCount};
 use crate::core::{
     models::{validate, DateLike, InvalidPaymentsError},
-    optimize::{brentq, newton_raphson},
+    optimize::{brentq, newton_raphson_2},
     utils::{self, fast_pow},
 };
 
@@ -15,16 +15,10 @@ pub fn xirr(
 
     let deltas = &day_count_factor(dates, day_count);
 
-    let f = |rate| {
-        if rate <= -1.0 {
-            // bound newton_raphson
-            return f64::INFINITY;
-        }
-        xnpv_result(amounts, deltas, rate)
-    };
-    let df = |rate| xnpv_result_deriv(amounts, deltas, rate);
+    let f = |rate| xnpv_result(amounts, deltas, rate);
+    let fd = |rate| xnpv_result_with_deriv(amounts, deltas, rate);
 
-    let rate = newton_raphson(guess.unwrap_or(0.1), &f, &df);
+    let rate = newton_raphson_2(guess.unwrap_or(0.1), &fd);
 
     if utils::is_a_good_rate(rate, f) {
         return Ok(rate);
@@ -39,7 +33,7 @@ pub fn xirr(
     let mut step = 0.01;
     let mut guess = -0.99999999999999;
     while guess < 1.0 {
-        let rate = newton_raphson(guess, &f, &df);
+        let rate = newton_raphson_2(guess, &fd);
         if utils::is_a_good_rate(rate, f) {
             return Ok(rate);
         }
@@ -86,15 +80,32 @@ fn day_count_factor(dates: &[DateLike], day_count: Option<DayCount>) -> Vec<f64>
 
 // \sum_{i=1}^n \frac{P_i}{(1 + rate)^{(d_i - d_0)/365}}
 fn xnpv_result(payments: &[f64], deltas: &[f64], rate: f64) -> f64 {
+    if rate <= -1.0 {
+        // bound newton_raphson
+        return f64::INFINITY;
+    }
     payments.iter().zip(deltas).map(|(p, &e)| p * fast_pow(1.0 + rate, -e)).sum()
 }
 
 // XNPV first derivative
 // \sum_{i=1}^n P_i * (d_0 - d_i) / 365 * (1 + rate)^{((d_0 - d_i)/365 - 1)}}
 // simplify in order to reuse cached deltas (d_i - d_0)/365
 // \sum_{i=1}^n \frac{P_i * -(d_i - d_0) / 365}{(1 + rate)^{((d_i - d_0)/365 + 1)}}
-fn xnpv_result_deriv(payments: &[f64], deltas: &[f64], rate: f64) -> f64 {
-    payments.iter().zip(deltas).map(|(p, e)| p * -e * fast_pow(1.0 + rate, -e - 1.0)).sum()
+// fn xnpv_result_deriv(payments: &[f64], deltas: &[f64], rate: f64) -> f64 {
+//     payments.iter().zip(deltas).map(|(p, e)| p * -e * fast_pow(1.0 + rate, -e - 1.0)).sum()
+// }
+
+fn xnpv_result_with_deriv(payments: &[f64], deltas: &[f64], rate: f64) -> (f64, f64) {
+    if rate <= -1.0 {
+        return (f64::INFINITY, f64::INFINITY);
+    }
+    // pow is an expensive function.
+    // we can re-use the result of pow for derivative calculation
+    payments.iter().zip(deltas).fold((0.0, 0.0), |acc, (p, e)| {
+        let y0 = p * fast_pow(1.0 + rate, -e);
+        let y1 = y0 * -e / (1.0 + rate);
+        (acc.0 + y0, acc.1 + y1)
+    })
 }
 
 #[cfg(test)]