arima_gpu.nu

arima_gpu.nu — the batch evaluator on a GPU, bit-for-bit the CPU's.

A fit asks for its likelihoods in rounds (src/arima.nu: the optimizer as a state machine), and K fits at once ask for K rounds together. This module answers such a round with one kernel launch: every item — a series, a method, a parameter vector — is one thread that runs the same filter the CPU runs, in the same order, with every operation rounded once (__dadd_rn and friends, never fused), so the value it returns is the value the CPU would have returned. The gain is the breadth: many series (a detector with a model per feature), or many candidate orders, in the time of one.

What stays on the host, and why: the transform from raw parameters to polynomials and the start of the model's stationary covariance (small, and they use libm — the device's log and exp are not the host's); and the sum of the innovation variances' logarithms, which the kernel hands back as the variances themselves. Only the filter's arithmetic — the Chandrasekhar recursions, O(n · r) — runs on the device.

Selection: any gpu backend — a CUDA device, or the package's host C++ backend (NURL_GPU=cpu) — gives the same numbers, so the choice is pace, never result. arima_fit_many_gpu falls back to the threaded CPU evaluator when no kit opens.

( arima_gpu_available ) → b ( arima_fit_many_gpu series spec method ) → ( Vec *ArimaModel ) ( arima_eval_gpu kit items ctxs out ) the evaluator itself

API

@ _ag_kernel_src → s

── the kernels ───────────────────────────────────────────────────────

meta, per item, 8 long longs: [0] w_off [1] n [2] r (ML) / p' (CSS) [3] f_off (ML) / ncond (CSS) [4] rmax (ML) / scratch_off (CSS) [5] poly_off [6] p0_off [7] q' (CSS)

: ~ b g_ag_trace F

Trace each round's phases (host preparation, the launches, the folds) to stderr — for tuning, off by default.

@ arima_gpu_set_trace b on → v

: ~ i g_ag_t_prep 0

: ~ i g_ag_t_run 0

: ~ i g_ag_t_post 0

: ~ i g_ag_rounds 0

@ arima_gpu_trace_report → v

@ arima_gpu_available → b

@ arima_eval_gpu * GpuKit kit ( Vec ArimaEvalItem ) items ( Vec ArimaCtx ) ctxs ( Vec f ) out → v

One round on the device: the ML items in one launch, the CSS items in another, each item's number folded on the host from what came back.

@ arima_fit_many_gpu ( Vec ( Vec f ) ) series ArimaSpec sp i method → ( Vec * ArimaModel )

K series fitted together with the device answering every round; the threaded CPU evaluator when no device (nor the CPU backend) opens.

---

arima.nu

arima — seasonal ARIMA forecasting for NURL: exact, fast, streaming.

A SARIMA(p, d, q)(P, D, Q)_s model in the form R's arima() and statsmodels' ARIMA estimate: the series is differenced d times and D times at the seasonal lag, the remainder is a stationary ARMA whose seasonal factors multiply the plain ones, and an optional mean sits under it all (only when nothing is differenced — a mean under a differenced series is a drift, which is a regressor, not a mean).

φ(B) Φ(B^s) (1 − B)^d (1 − B^s)^D (y_t − μ) = θ(B) Θ(B^s) ε_t

Estimation is CSS-ML: conditional sum of squares first (cheap, and a good place to start), then exact maximum likelihood by the Kalman filter over the state-space form (Harvey), the initial state covariance being the model's own stationary covariance — solved exactly, by the doubling recursion — so the likelihood is the true Gaussian likelihood of the differenced series, not an approximation that depends on how the series began. σ² is concentrated out. The optimizer is BFGS with a backtracking line search over parameters transformed (Jones 1980, the PACF map R uses) so every AR polynomial it tries is stationary and every MA polynomial invertible.

Forecasting and streaming run on the FULL state-space model — the ARMA state with the differencing folded in, the way R's makeARIMA builds it — filtered over the raw series with a diffuse start for the differencing states. So a forecast's mean and its standard error come from the model's own state at the end of the data, and one new observation is one Kalman step: arima_update absorbs it in O(r²), reports the innovation it made and how surprising it was, and the next forecast starts from there. That is the streaming shape the anomaly service wants: a model trained once, kept current point by point, and refitted on a schedule.

Surface (see README.md for the story): ( arima_spec p d q ) / ( arima_spec_seasonal p d q P D Q s ) → ArimaSpec ( arima_fit y spec ) → ArimaModel CSS-ML, the default ( arima_fit_method y spec m ) → ArimaModel m = ARIMA_CSS | ARIMA_ML ( arima_auto y s ) → *ArimaModel stepwise order search by AICc ( arima_forecast m h ) → ArimaForecast h means and standard errors ( arima_update m y ) → ArimaStep one observation in: innovation, variance, z ( arima_coef m ) → Json coefficients, σ², log-likelihood, AIC… ( arima_to_json m ) / ( arima_from_json s ) persistence, bit-exact ( arima_free m )

Pure NURL, no dependencies beyond the stdlib; src/arima_gpu.nu adds batched fitting on a GPU through the gpu package for the case of many series or many candidate orders at once.

API

: i ARIMA_CSS 0

: i ARIMA_ML 1

: f ARIMA_KAPPA 1000000.0

The diffuse prior on a differencing state: R's kappa.

: i ARIMA_MAX_ITER 200

Where the optimizer stops.

: f ARIMA_TOL 0.00000001

: f ARIMA_LOG_2PI 1.8378770664093453

: i ARIMA_PAR_WORK 200000

The numerical gradient's evaluations run on threads once one evaluation is worth a thread: this many multiply-adds per likelihood.

: ArimaSpec

: ArimaSpec {
    i p
    i d
    i q
    i P
    i D
    i Q
    i s  // seasonal period; 0 or 1 = none
    b mean  // fit a mean (only honoured when d + D = 0)
}

@ arima_spec i p i d i q → ArimaSpec

@ arima_spec_seasonal i p i d i q i P i D i Q i s → ArimaSpec

@ arima_spec_with_mean ArimaSpec sp b mean → ArimaSpec

@ _ar_expand_ar ( Vec f ) phi ( Vec f ) sphi i s → ( Vec f )

The AR side as "1 − Σ φ B^k" coefficients: φ(B)Φ(B^s) expanded, the signs such that the result is the φ_k of the expanded polynomial.

@ _ar_expand_ma ( Vec f ) theta ( Vec f ) stheta i s → ( Vec f )

The MA side: θ(B)Θ(B^s) expanded, "1 + Σ θ B^k" convention as given.

@ _ar_delta i d i D i s → ( Vec f )

The differencing polynomial (1 − B)^d (1 − B^s)^D as "1 − Σ δk B^k": returns δ (length d + s·D), so that yt = w_t + Σ δk y{t−k}.

@ arima_difference ( Vec f ) y i d i D i s → ( Vec f )

Δ^d Δ_s^D y: the first d + s·D values are consumed.

@ _ar_partrans ( Vec f ) raw i off i n → ( Vec f )

@ _ar_invpartrans ( Vec f ) phi ( Vec f ) raw i off → b

Inverse: AR coefficients → raw. Returns F when the polynomial is not stationary (a partial autocorrelation reaches 1 in magnitude).

: ArimaCoef

: ArimaCoef {
    ( Vec f ) phi
    ( Vec f ) theta
    ( Vec f ) sphi
    ( Vec f ) stheta
    f mu
}

@ _ar_coef_new ArimaSpec sp → ArimaCoef

@ _ar_coef_free ArimaCoef c → v

@ _ar_coef_of_raw ArimaSpec sp ( Vec f ) raw → ArimaCoef

Raw (transformed) vector → coefficients. Layout: φ, θ, Φ, Θ, μ.

@ _ar_at ( Vec f ) v i idx → f

: ArimaSS

: ArimaSS {
    i r
    i nd
    i rd
    ( Vec f ) phi  // r, the expanded AR coefficients padded with zeros
    ( Vec f ) theta  // r, θ_0 = 1 then the expanded MA padded
    ( Vec f ) delta  // nd
    ( Vec f ) a  // rd
    ( Vec f ) pm  // rd × rd, row-major
    ( Vec f ) scratch  // rd × rd
    ( Vec f ) scratch2  // rd × rd
    ( Vec f ) pz  // rd
    ( Vec f ) prev  // rd × rd, the covariance a step ago (the steady-state test)
    ( Vec f ) kg  // rd, the gain once the covariance has converged
    ( Vec f ) fz  // 2: [the innovation variance at convergence, 1.0 once converged]
}

@ _ar_ss_steady ArimaSS ss → b

Has the covariance recursion converged (see arstep)?

@ _ar_ss_free ArimaSS ss → v

@ _ar_ss_new ( Vec f ) ar ( Vec f ) ma ( Vec f ) delta → ArimaSS

Build the form from expanded polynomials (ar: "1 − Σ φ B^k" φ's; ma: "1 + Σ θ B^k" θ's) and the differencing δ.

: ArimaStep

: ArimaStep {
    f innovation  // y − its one-step forecast
    f variance  // the forecast's variance, σ² units applied by the caller
    f predicted  // the one-step forecast that was made
}

@ _ar_step ArimaSS ss f y → ArimaStep

One filter step: observe y, update, predict the next. Returns the innovation and its variance in σ² = 1 units; F ≤ 0 makes the step report a negative variance (the caller treats that as failure). One filter step of the full model. The covariance recursion converges (the model is time-invariant); once the step moved it by no more than 10⁻¹⁴ (1 + F) in any element the gain is fixed — the steady state — and a step is O(r_d): the innovation, the state moved by the gain, the transition. Before that, the O(r_d²) covariance form. The two agree to the last bit with what the covariance form would go on producing, short of the increments it stopped adding.

@ _ar_solve ( Vec f ) M ( Vec f ) b i n → b

Solve M x = b for a small dense system by Gaussian elimination with partial pivoting; M and b are overwritten, x lands in b. F when singular.

@ _ar_autocov ( Vec f ) ar ( Vec f ) ma i m ( Vec f ) gamma ( Vec f ) psi → b

The autocovariances γ(0..m−1) of the ARMA(p', q') with expanded polynomials (σ² = 1): the Yule–Walker system for γ(0..p'), then the recursion. Also the ψ weights ψ(0..m−1). F when the system is singular.

@ _ar_init_cov ArimaSS ss → b

The stationary covariance of the ARMA block, from the autocovariances: with a_t[i] = Σm (φ{i+1+m} y_{t−1−m} + θ{i+m} ε{t−m}), P = Φ Γ Φᵀ + Φ C Θᵀ + Θ Cᵀ Φᵀ + Θ Θᵀ, Φ[i][m] = φ{i+1+m}, Θ[i][m] = θ{i+m}, Γ[m][l] = γ(|m−l|), C[m][l] = ψ_{l−1−m} (l > m). Three r³ products, no iteration — the doubling recursion needed dozens of them for a seasonal polynomial's roots close to the circle. Writes the ARMA block of pm; the differencing block gets the diffuse prior.

@ _ar_init_cov_doubling ArimaSS ss → b

The same by doubling (P_{k+1} = P_k + A_k P_k A_kᵀ, A_{k+1} = A_k²): kept as the independent check of the closed form.

: ArimaLik

: ArimaLik {
    b ok
    f loglik
    f sigma2
    i n_used
}

: ArimaArma

: ArimaArma {
    i r
    ( Vec f ) phi
    ( Vec f ) theta
}

The ARMA block alone, padded to the state's width: φ1..φr and θ0 = 1, θ1..θ_{r−1} — what the exact likelihood works from.

@ _ar_arma_new ( Vec f ) ar ( Vec f ) ma → ArimaArma

@ _ar_arma_free ArimaArma a → v

@ _ar_init_col i r ( Vec f ) phiv ( Vec f ) thv ( Vec f ) out → b

The first column of the stationary covariance, P e₀, in O(r²): the four terms of arinit_cov applied to the unit vector instead of multiplied out — Φ(Γ Φᵀe₀) + Φ(C Θᵀe₀) + Θ(Bᵀe₀) + Θ(Θᵀe₀) with B = ΦC. It is all the Chandrasekhar recursion needs of P.

@ _ar_filter_arma ( Vec f ) phiv ( Vec f ) thv ( Vec f ) col ( Vec f ) w f mu → ArimaLik

The exact Gaussian likelihood of the stationary series w (mean μ) under the ARMA block, by the Chandrasekhar recursions (Morf, Sidhu & Kailath 1974; Herbst 2015 for this form). The filter's covariance recursion, started at the stationary P, moves by a rank-one increment P_{t+1} − P_t = W_t M_t W_tᵀ, and that increment has its own recursion — F_{t+1} = F_t + (Z W_t)² M_t, K_{t+1} = (K_t F_t + T W_t M_t Z W_t)/F_{t+1}, W_{t+1} = (T − K_{t+1} Z) W_t, M_{t+1} = M_t + (M_t Z W_t)²/F_t — so a step costs O(r) where the covariance form costs O(r²): the same innovations and variances, arrived at without P. col is P e₀ (the start needs nothing else: F₁ = P₀₀, K₁ = T P e₀ / F₁, W₁ = K₁, M₁ = −F₁). Once the increment is below 10⁻¹⁴ (1 + F) in every element the gain is fixed — the steady state — and only the state moves.

@ _ar_lik_ml_from f ssq f sumlog i n → ArimaLik

The exact log-likelihood from the filter's two sums — one formula for every route that produces them.

@ _ar_lik_css_from f ssq i nu → ArimaLik

The conditional-sum-of-squares "log-likelihood" from its sum.

: ArimaObj

: ArimaObj {
    ArimaSpec sp
    ( Vec f ) w
    i method
    i ncond
    i evals
}

What the objective needs to see: the spec, the differenced series, the method (and CSS's conditioning count), and a count of evaluations.

: ArimaJob

: ArimaJob {
    ArimaSpec sp
    ( Vec f ) w
    i method
    i ncond
    ( Vec f ) raw
    ( Vec f ) out
    i idx
    i kind  // 0 raw objective, 1 natural objective, 2/3 prepare (raw/natural, ML), 4/5 prepare (raw/natural, CSS)
    i extra  // a prepare job's *ArimaPrep
}

One evaluation as a job: its parameters (raw, or — kind 1 — natural coefficients for the Hessian), its slot in the answers.

: ArimaLane

: ArimaLane {
    ( Vec i ) jobs
    i lane
    i stride
}

A worker: every job whose index ≡ lane (mod stride).

@ _ar_jobs_run ( Vec i ) jobs b par → v

Run every job: on a pool of __ar_threads workers striding the list when par, else in place. Frees the jobs.

@ _ar_jobs_free ( Vec i ) jobs → v

Free the jobs of a batch (their results have been read).

: ArimaPrep

: ArimaPrep {
    b ok
    i r
    ( Vec f ) phi
    ( Vec f ) theta
    ( Vec f ) p0  // ML: P e₀, r values
    ( Vec f ) ar
    ( Vec f ) ma
    f mu
}

A state-space form at a parameter point, prepared for a device: the expanded polynomials, the padded state vectors and, for ML, the first column of the stationary covariance. What the kernel cannot compute itself.

@ _ar_prep_free * ArimaPrep p → v

@ _ar_prep_at ArimaSpec sp ( Vec f ) raw b natural b with_cov → *ArimaPrep

Prepare at raw (transformed parameters when natural = F, natural coefficients when T); with_cov adds the stationary covariance.

@ _ar_job_new ArimaSpec sp ( Vec f ) w i method i ncond ( Vec f ) raw ( Vec f ) out i idx i kind → i

: i ARIMA_PH_START 0

: i ARIMA_PH_GRAD 1

: i ARIMA_PH_TRIAL 2

: ArimaBfgs

: ArimaBfgs {
    i k
    ( Vec f ) raw
    ( Vec f ) g
    ( Vec f ) d
    ( Vec f ) s
    ( Vec f ) yv
    ( Vec f ) hy
    ( Vec f ) trial
    ( Vec f ) H
    ( Vec f ) hs  // the finite-difference steps of the pending stencil
    f f0
    f fn
    f alpha
    f gd
    i tries
    i iter
    i phase
    b converged
    b done
}

: ArimaOpt

: ArimaOpt {
    b converged
    f value
    i iterations
}

: ArimaModel

: ArimaModel {
    ArimaSpec spec
    ArimaCoef coef
    f sigma2
    f loglik
    f aic
    f aicc
    f bic
    i n  // observations the model has absorbed (fit + updates)
    i n_fit  // observations the coefficients were fitted on
    i n_used  // differenced observations in the likelihood
    i method
    b converged
    i iterations
    i evals
    ArimaSS ss  // the full model's state after the last observation
    ( Vec f ) se  // standard errors, coefficient order (may be NaN)
    f last_innovation
    f last_variance
    f last_predicted
    ( Vec i ) xper  // Fourier periods in rows (empty: no regressors) — see arima_fit_harmonic
    i xk  // harmonics per period
    ( Vec f ) xcoef  // 1 + 2·xk·|xper|: the intercept, then per period and harmonic the sine and cosine weights
    i xt  // rows since the fit's origin: the index the next observation has in the regressors
}

@ arima_free * ArimaModel m → v

@ arima_spec_of * ArimaModel m → ArimaSpec

@ arima_sigma2 * ArimaModel m → f

@ arima_loglik * ArimaModel m → f

@ arima_aic * ArimaModel m → f

@ arima_aicc * ArimaModel m → f

@ arima_n * ArimaModel m → i

@ arima_converged * ArimaModel m → b

@ arima_phi * ArimaModel m → ( Vec f )

@ arima_theta * ArimaModel m → ( Vec f )

@ arima_sphi * ArimaModel m → ( Vec f )

@ arima_stheta * ArimaModel m → ( Vec f )

@ arima_mu * ArimaModel m → f

@ _ar_coef_of_natural ArimaSpec sp ( Vec f ) x → ArimaCoef

Natural coefficients (φ, θ, Φ, Θ, μ in one vector) → the bundle.

@ _ar_natural_of_coef ArimaSpec sp ArimaCoef mc → ( Vec f )

The bundle → the natural vector.

@ arima_fit_method ( Vec f ) y ArimaSpec sp0 i method → *ArimaModel

: ArimaEvalItem

: ArimaEvalItem {
    i ctx
    ( Vec f ) raw
    i kind  // 0 = transformed parameters, 1 = natural coefficients
}

One evaluation of a batch: which context (series + method), at which parameters.

: ArimaCtx

: ArimaCtx {
    ArimaSpec sp
    ( Vec f ) w
    i method
    i ncond
}

A context: the differenced series, the specification, the method (and CSS's conditioning count, 0 = the model's own order).

@ arima_eval_cpu ( Vec ArimaEvalItem ) items ( Vec ArimaCtx ) ctxs ( Vec f ) out → v

The threaded CPU evaluator: every item on a thread of its own, __ar_threads at a time, when one is worth it.

: ArimaFitState

: ArimaFitState {
    * ArimaBfgs st
    i stage  // 0 = CSS running, 1 = ML running, 2 = finished
    i iters
    i evals
    b converged
    i req_at  // where this model's requests start in the current batch
    i req_n
}

@ arima_fit_many_with ( Vec ( Vec f ) ) series ArimaSpec sp0 i method ( @ v ( Vec ArimaEvalItem ) ( Vec ArimaCtx ) ( Vec f ) ) evaluator → ( Vec * ArimaModel )

Fit series (each a raw series) under sp by method, the evaluations batched through evaluator.

@ _ar_geti ( Vec i ) v i k → i

@ arima_fit_many ( Vec ( Vec f ) ) series ArimaSpec sp i method → ( Vec * ArimaModel )

K series fitted together on the CPU's threads.

@ arima_models_free ( Vec * ArimaModel ) ms → v

@ arima_fit ( Vec f ) y ArimaSpec sp → *ArimaModel

: ArimaForecast

: ArimaForecast {
    ( Vec f ) mean
    ( Vec f ) se
}

@ arima_forecast_free ArimaForecast fc → v

@ arima_forecast * ArimaModel m i h → ArimaForecast

h steps ahead from the model's current state: means and standard errors (σ² applied). The state is left where it was.

: ArimaUpdate

: ArimaUpdate {
    f predicted  // what the model expected
    f innovation  // y − predicted
    f variance  // the forecast variance (σ² applied)
    f z  // innovation / √variance: how surprising y was
}

@ arima_update * ArimaModel m f y → ArimaUpdate

One new observation: a Kalman step on the full model. The coefficients do not change; refit when the schedule says so.

A NaN y is a missing observation: the step is the filter's time update alone — the state moves on, its uncertainty grows, nothing is learned — so a gap in a stream costs a tick of the clock and not a restart. The answer then carries what the model predicted and the variance it would have judged an observation by; innovation and z are NaN, because there was nothing to be surprised by.

@ arima_restart * ArimaModel m → v

Forget the observations: the state goes back to where a freshly fitted model's stands before its first point (the stationary covariance for the ARMA part, diffuse for the differencing), n to 0. The coefficients stay. Feed the series again with arima_update and the state after n points is the state a fit over them would have left — which is how a caller replays a stored history through a model whose state has moved past it.

@ arima_restart_at * ArimaModel m i t0 → v

The same, with the regressors' clock set: t0 is the row the next observation has, counted from the fit's origin (negative for rows before it) — a replay that begins elsewhere than the fit did keeps the seasonal's phase. A model without regressors ignores it.

@ arima_clone * ArimaModel m → *ArimaModel

A deep copy: coefficients, fit statistics, standard errors and the state, so the copy can be stepped without moving the original.

@ arima_fit_harmonic ( Vec f ) y ( Vec i ) periods i k ArimaSpec sp i method → *ArimaModel

Fit y as Fourier terms of periods (rows) with k harmonics each, by least squares, then the ARIMA sp on the residuals by method, and attach the terms to the model: its updates and forecasts carry them, its state and statistics are the residual model's. t = 0 is the first row of y.

@ arima_auto_harmonic ( Vec f ) y ( Vec i ) periods i k i s → *ArimaModel

The same with the residual model's order chosen by the stepwise search (s its season, 0 for none — the periods carry the long ones).

@ arima_kpss ( Vec f ) x → f

The KPSS statistic for level stationarity of x, with the Bartlett long-run variance over l = 4 (n/100)^{1/4} lags (Kwiatkowski et al.).

@ arima_ndiffs ( Vec f ) y → i

Differences needed for level stationarity: KPSS at 5 % (0.463), at most two.

@ arima_acf ( Vec f ) x i k → f

The sample autocorrelation of x at lag k.

@ arima_nsdiffs ( Vec f ) y i s i d → i

A seasonal difference when the season carries the series: the autocorrelation at lag s of the (first-differenced, if d > 0) series above 0.5 — the rule of thumb behind the older auto.arima.

@ arima_auto ( Vec f ) y i s → *ArimaModel

Stepwise search; s is the season (0 = none), d and D chosen by the tests above. Returns the best model found.

: i ARIMA_SCREEN_N 150

Candidates are screened by conditional sum of squares — the way auto.arima approximates — whenever the series is longer than ARIMA_SCREEN_N points or the season longer than ARIMA_SCREEN_S (R's rule: approximation = n > 150 | frequency > 12), and the winner is then refitted by ML. Below that, every candidate gets the exact likelihood. The rule is about the search, not one fit: a daily season on hourly data makes a candidate's state 50 wide, and the stepwise search evaluates thousands of likelihoods — 42 s exactly against 0.3 s screened, for the same chosen order.

: i ARIMA_SCREEN_S 12

@ arima_auto_d ( Vec f ) y i s i d i D → *ArimaModel

@ _ar_jarr ( Vec f ) v → Json

@ arima_coef * ArimaModel m → Json

Coefficients and fit statistics as JSON, for a report or a table.

@ arima_to_json * ArimaModel m → String

@ arima_from_json s src → ?*ArimaModel

---

main.nu

arima — the CLI: fit or select a model on a column of numbers, forecast.

arima fit FILE --order p,d,q [--seasonal P,D,Q,s] [--mean] [--css] [--horizon h] arima auto FILE [--season s] [--horizon h]

FILE holds one number per line (a CSV's first column is taken; a header line that is not a number is skipped). The answer is JSON: the coefficients and fit statistics (arima_coef) and, with --horizon, the forecast means and standard errors.

API

@ usage → i

@ read_series s path → ( Vec f )

Read the first numeric column of a text file.

@ ints_of String spec ( Vec i ) out → v

"1,1,1" → up to 4 integers (missing ones 0).

@ arg_after ( Vec String ) args s flag → String

@ has_flag ( Vec String ) args s flag → b

@ report * ArimaModel m i h → v

@ main → i

@ vec_get_i ( Vec i ) v i k → i