perf(linreg): incremental O(1) OLS for LinearRegression and LinRegSlope (R2)

`LinearRegression::fit` and `LinRegSlope::update` previously iterated the
full `period`-window on every tick to recompute `Σy` and `Σxy` from
scratch — O(period) per update, in violation of the `Indicator` trait's
O(1) contract. `LinRegAngle` inherits the cost transitively because it
delegates to `LinRegSlope`.

This commit slides the OLS state in closed form. The constant terms
(`Σx`, `Σxx`, the denominator `n·Σxx − (Σx)²`) were already precomputed
in `new`. The new running state is:

- `sum_y: f64` — running sum of the values currently in the window.
- `sum_xy: f64` — running Σ(x · y) where `x` is the position of each
  value inside the trailing window (`0` for the oldest, `n−1` for the
  newest).

On every push, when the window is already full the front value `y₀` is
popped and the indices of every remaining value shift down by 1; the
identity

    new_Σxy = old_Σxy − old_Σy + y₀

closes the slide in O(1). The new value is then pushed at position `k`
(the current length before the push), contributing `k · new_value` to
`sum_xy` and `new_value` to `sum_y`. The output is the same TA-Lib OLS
formula evaluated against the incremental accumulators.

Behaviour is unchanged: same per-tick values, same warmup, same NaN
semantics. Two new tests compare the O(1) result bar-by-bar against a
fresh O(n) refit on a noisy ramp (sliding-phase dominated), a step
function (large pop/push deltas), and constants (tests floating-point
drift) — agreement is within `1e-9`.

`LinRegAngle` benefits automatically through its `LinRegSlope` field.
This commit is contained in:
kingchenc
2026-05-23 10:36:45 +02:00
parent b003321562
commit 2aef8c8db5
3 changed files with 192 additions and 28 deletions
+11
View File
@@ -24,6 +24,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
site makes the invariant explicit.
### Changed
- `LinearRegression`, `LinRegSlope` and `LinRegAngle` (via composition over
`LinRegSlope`) now run their rolling ordinary-least-squares fit
**incrementally** in O(1) per update (audit finding R2). Previously every
tick refit the line from scratch in O(period). The OLS denominators (`Σx`
and `Σxx`) depend only on `period`, so they were already precomputed; this
release adds running `Σy` and `Σxy` accumulators and slides them in closed
form via the identity
`new_Σxy = old_Σxy old_Σy + popped_y₀` (then `Σxy += (n 1) · new_value`
and `Σy += new_value`). New per-bar equivalence tests compare the O(1)
output against a fresh O(n) refit on noisy ramps, step functions, and
constants — values agree to within 1e-9.
- Fuzz suite expanded from 2 indicators to the full catalogue (audit finding
R9). The existing `indicator_update` target now exercises every scalar-input
indicator (~33 classes including MACD and Bollinger Bands); a new
+105 -18
View File
@@ -19,8 +19,20 @@ use crate::traits::Indicator;
///
/// This is TA-Lib's `LINEARREG`: a smoothed price that lags less than an SMA
/// because it extrapolates the *local trend* forward to the current bar
/// instead of averaging it away. The `Σx` terms depend only on `period`, so
/// they are computed once; each `update` is O(period).
/// instead of averaging it away.
///
/// Each `update` is O(1): the `Σx` and `Σxx` terms depend only on `period` and
/// are precomputed once, while `Σy` and `Σxy` are maintained incrementally as
/// the window slides. The closed-form sliding-window identity for
/// `x = 0, 1, …, period 1` is
///
/// ```text
/// new_sum_xy = old_sum_xy old_sum_y + popped_y0 // index shift by 1
/// new_sum_y = old_sum_y popped_y0
/// // then push the new value at index n1:
/// sum_xy += (n 1) · new_value
/// sum_y += new_value
/// ```
///
/// # Example
///
@@ -38,8 +50,16 @@ use crate::traits::Indicator;
pub struct LinearRegression {
period: usize,
window: VecDeque<f64>,
/// Closed form of `Σx` over `x = 0, 1, …, period 1` — constant in `period`.
sum_x: f64,
/// Closed form of `n · Σxx (Σx)²` — constant in `period`, the OLS
/// denominator.
denom: f64,
/// Running sum of the values currently in the window.
sum_y: f64,
/// Running `Σ(x · y)` where `x` is the position of each value within the
/// trailing window (`0` for the oldest, `period 1` for the newest).
sum_xy: f64,
}
impl LinearRegression {
@@ -63,6 +83,8 @@ impl LinearRegression {
window: VecDeque::with_capacity(period),
sum_x,
denom: n * sum_xx - sum_x * sum_x,
sum_y: 0.0,
sum_xy: 0.0,
})
}
@@ -70,20 +92,6 @@ impl LinearRegression {
pub const fn period(&self) -> usize {
self.period
}
/// Ordinary-least-squares `(slope, endpoint)` over the current full window.
fn fit(&self) -> (f64, f64) {
let n = self.period as f64;
let mut sum_y = 0.0;
let mut sum_xy = 0.0;
for (x, &y) in self.window.iter().enumerate() {
sum_y += y;
sum_xy += x as f64 * y;
}
let slope = (n * sum_xy - self.sum_x * sum_y) / self.denom;
let intercept = (sum_y - slope * self.sum_x) / n;
(slope, intercept + slope * (n - 1.0))
}
}
impl Indicator for LinearRegression {
@@ -92,17 +100,35 @@ impl Indicator for LinearRegression {
fn update(&mut self, value: f64) -> Option<f64> {
if self.window.len() == self.period {
self.window.pop_front();
// Sliding phase: pop the oldest, then shift every remaining index
// down by 1 in the running `sum_xy`. The identity
// Σ((i 1) · y_i for i = 1..n1) = Σ(i · y_i) Σ(y_i) + y_0
// gives the closed-form update below.
let y0 = self.window.pop_front().expect("non-empty");
self.sum_xy = self.sum_xy - self.sum_y + y0;
self.sum_y -= y0;
}
// Append at position `k = current length` before the push. During
// warmup `k` ranges over `0..period 1`; once the window is full it
// is always `period 1`.
let k = self.window.len() as f64;
self.window.push_back(value);
self.sum_y += value;
self.sum_xy += k * value;
if self.window.len() < self.period {
return None;
}
Some(self.fit().1)
let n = self.period as f64;
let slope = (n * self.sum_xy - self.sum_x * self.sum_y) / self.denom;
let intercept = (self.sum_y - slope * self.sum_x) / n;
Some(intercept + slope * (n - 1.0))
}
fn reset(&mut self) {
self.window.clear();
self.sum_y = 0.0;
self.sum_xy = 0.0;
}
fn warmup_period(&self) -> usize {
@@ -195,4 +221,65 @@ mod tests {
prices.iter().map(|x| b.update(*x)).collect::<Vec<_>>()
);
}
/// Incremental OLS equivalence: the O(1) implementation must agree to
/// `1e-9` with a fresh-from-scratch O(n) refit on every bar, on inputs
/// chosen to stress every code path: a noisy ramp (sliding phase
/// dominates), a step function (the new value differs sharply from the
/// popped one), and constants (the floating-point accumulators must not
/// drift).
#[test]
fn incremental_matches_naive_fit_bar_by_bar() {
fn naive_endpoint(window: &[f64]) -> f64 {
let n = window.len() as f64;
let mut sum_y = 0.0;
let mut sum_xy = 0.0;
let mut sum_x = 0.0;
let mut sum_xx = 0.0;
for (i, &y) in window.iter().enumerate() {
let x = i as f64;
sum_y += y;
sum_xy += x * y;
sum_x += x;
sum_xx += x * x;
}
let denom = n * sum_xx - sum_x * sum_x;
let slope = (n * sum_xy - sum_x * sum_y) / denom;
let intercept = (sum_y - slope * sum_x) / n;
intercept + slope * (n - 1.0)
}
fn check(prices: &[f64], period: usize) {
let mut lr = LinearRegression::new(period).unwrap();
for (t, p) in prices.iter().enumerate() {
let streaming = lr.update(*p);
if t + 1 >= period {
let lo = t + 1 - period;
let expected = naive_endpoint(&prices[lo..=t]);
let got = streaming.expect("warmed up");
assert!(
(got - expected).abs() < 1e-9,
"endpoint diverges at t={t}, period={period}: got={got}, expected={expected}",
);
}
}
}
let noisy_ramp: Vec<f64> = (0..120)
.map(|i| 100.0 + f64::from(i) * 0.5 + (f64::from(i) * 0.7).sin() * 3.0)
.collect();
check(&noisy_ramp, 5);
check(&noisy_ramp, 14);
check(&noisy_ramp, 30);
let mut step = vec![1.0; 30];
step.extend(std::iter::repeat_n(100.0, 30));
step.extend(std::iter::repeat_n(0.001, 30));
check(&step, 5);
check(&step, 14);
let constant = vec![42.0; 50];
check(&constant, 8);
check(&constant, 25);
}
}
@@ -17,8 +17,12 @@ use crate::traits::Indicator;
/// This is TA-Lib's `LINEARREG_SLOPE`: a momentum-like reading of how steeply
/// price is trending over the window — positive while it rises, negative
/// while it falls, near zero when it is flat — without the band-pass quirks
/// of a difference-based oscillator. The `Σx` terms depend only on `period`,
/// so they are computed once; each `update` is O(period).
/// of a difference-based oscillator.
///
/// Each `update` is O(1): the same incremental OLS state as
/// [`LinearRegression`](crate::LinearRegression) is maintained — `Σx` and
/// `Σxx` are precomputed once from `period`, while `Σy` and `Σxy` are slid
/// forward in closed form on every push.
///
/// # Example
///
@@ -36,8 +40,14 @@ use crate::traits::Indicator;
pub struct LinRegSlope {
period: usize,
window: VecDeque<f64>,
/// Closed form of `Σx` over `x = 0, 1, …, period 1` — constant in `period`.
sum_x: f64,
/// Closed form of `n · Σxx (Σx)²` — constant in `period`.
denom: f64,
/// Running sum of the values currently in the window.
sum_y: f64,
/// Running `Σ(x · y)` where `x` is the position within the trailing window.
sum_xy: f64,
}
impl LinRegSlope {
@@ -61,6 +71,8 @@ impl LinRegSlope {
window: VecDeque::with_capacity(period),
sum_x,
denom: n * sum_xx - sum_x * sum_x,
sum_y: 0.0,
sum_xy: 0.0,
})
}
@@ -76,24 +88,30 @@ impl Indicator for LinRegSlope {
fn update(&mut self, value: f64) -> Option<f64> {
if self.window.len() == self.period {
self.window.pop_front();
// Sliding-window identity: when the window slides one step forward
// the indices `x` for every kept entry shift down by 1, so
// new_sum_xy = old_sum_xy old_sum_y + y0
// (`y0` is the popped front value).
let y0 = self.window.pop_front().expect("non-empty");
self.sum_xy = self.sum_xy - self.sum_y + y0;
self.sum_y -= y0;
}
let k = self.window.len() as f64;
self.window.push_back(value);
self.sum_y += value;
self.sum_xy += k * value;
if self.window.len() < self.period {
return None;
}
let n = self.period as f64;
let mut sum_y = 0.0;
let mut sum_xy = 0.0;
for (x, &y) in self.window.iter().enumerate() {
sum_y += y;
sum_xy += x as f64 * y;
}
Some((n * sum_xy - self.sum_x * sum_y) / self.denom)
Some((n * self.sum_xy - self.sum_x * self.sum_y) / self.denom)
}
fn reset(&mut self) {
self.window.clear();
self.sum_y = 0.0;
self.sum_xy = 0.0;
}
fn warmup_period(&self) -> usize {
@@ -192,4 +210,52 @@ mod tests {
prices.iter().map(|x| b.update(*x)).collect::<Vec<_>>()
);
}
/// Incremental OLS equivalence for the slope: the O(1) implementation must
/// agree bar-by-bar with a fresh-from-scratch O(n) refit, on a noisy ramp
/// (sliding-phase dominated) and a step function (large pop/push deltas).
#[test]
fn incremental_matches_naive_slope_bar_by_bar() {
fn naive_slope(window: &[f64]) -> f64 {
let n = window.len() as f64;
let mut sum_y = 0.0;
let mut sum_xy = 0.0;
let mut sum_x = 0.0;
let mut sum_xx = 0.0;
for (i, &y) in window.iter().enumerate() {
let x = i as f64;
sum_y += y;
sum_xy += x * y;
sum_x += x;
sum_xx += x * x;
}
(n * sum_xy - sum_x * sum_y) / (n * sum_xx - sum_x * sum_x)
}
fn check(prices: &[f64], period: usize) {
let mut ls = LinRegSlope::new(period).unwrap();
for (t, p) in prices.iter().enumerate() {
let streaming = ls.update(*p);
if t + 1 >= period {
let lo = t + 1 - period;
let expected = naive_slope(&prices[lo..=t]);
let got = streaming.expect("warmed up");
assert!(
(got - expected).abs() < 1e-9,
"slope diverges at t={t}, period={period}: got={got}, expected={expected}",
);
}
}
}
let noisy_ramp: Vec<f64> = (0..120)
.map(|i| 100.0 + f64::from(i) * 0.5 + (f64::from(i) * 0.7).sin() * 3.0)
.collect();
check(&noisy_ramp, 5);
check(&noisy_ramp, 14);
let mut step = vec![1.0; 30];
step.extend(std::iter::repeat_n(100.0, 30));
check(&step, 7);
}
}