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 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 part that is 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 [4] scratch_off [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
}

@ _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_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
}

@ _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 evals
}

What the objective needs to see: the spec, the differenced series, the method, and a count of evaluations.

: ArimaJob

: ArimaJob {
    ArimaSpec sp
    ( Vec f ) w
    i method
    ( 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).

: ArimaPrep

: ArimaPrep {
    b ok
    i r
    ( Vec f ) phi
    ( Vec f ) theta
    ( Vec f ) p0
    ( 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 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.

: 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
}

@ 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
}

A context: the differenced series, the specification, the method.

@ 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.

@ 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_WORK 20000000

Candidates are screened by CSS when an exact likelihood of the largest of them would be expensive (n · r² beyond ARIMA_SCREEN_WORK — a weekly season on hourly data), the way auto.arima approximates; the winner is then refitted by ML.

@ 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