CPU Parallelism¶

pyscarcopula provides two complementary CPU parallelism levels:

native C++ threads inside one multivariate fit or numerical kernel;
worker processes for independent fits and rolling windows.

Parallelism is opt-in. If no thread argument is supplied, every multivariate model and numerical method uses exactly one native thread, independently of the process environment.

Enabling native threads¶

Pass a NumericalConfig to a fit when the complete fit should use native threads:

import numpy as np

from pyscarcopula import (
    EquicorrGaussianCopula,
    NumericalConfig,
    StochasticStudentCopula,
)

cfg = NumericalConfig(n_threads=4)

equicorr = EquicorrGaussianCopula(d=u.shape[1])
equicorr_result = equicorr.fit(
    u,
    method="scar-tm-ou",
    config=cfg,
)

student = StochasticStudentCopula(d=u.shape[1], R=R)
student_result = student.fit(
    u,
    method="mle",
    config=cfg,
)

assert equicorr_result.diagnostics["n_threads"] == 4
assert student_result.diagnostics["n_threads"] == 4

Methods that expose n_threads directly can opt in per call:

rows = equicorr.log_pdf_rows(u, r=0.25, n_threads=4)

conditional = student.sample_conditional(
    20_000,
    r=6.0,
    given={0: 0.25, 2: 0.8},
    rng=np.random.default_rng(2026),
    n_threads=4,
)

For native method arguments and NumericalConfig, valid explicit values are integers from 1 through 256. Boolean values, 0, None, non-integral values, and values above 256 are rejected. Process helpers use their documented None default only as an omitted-value sentinel and resolve it to the strict one-thread default.

Absolute one-thread default¶

The default is intentionally stricter than an automatic CPU-count policy:

from pyscarcopula import NumericalConfig

assert NumericalConfig().n_threads == 1

No environment variable can change this result. In particular, PYSCARCOPULA_NUM_THREADS is not supported or read. Native parallelism must be requested through an explicit method argument or NumericalConfig(n_threads=N).

With n_threads=1, the native thread pool is not created or consulted. This keeps sequential execution safe for an outer rolling-window executor and avoids hidden background workers.

Runtime work and the resident pool¶

The native runtime uses the following vocabulary when it plans one call:

Symbol	Meaning
`R`	Numerical partials required by the algorithm, such as `min(N, 64)` for a Factor Student joint reduction
`B`	Stable logical portions, including their ranges and IDs; some reductions use `B = R`
`W`	Executor budget for this call after validation and any kernel-specific cost limit, never above the explicit thread request
`J`	Queued runners for a non-empty prepared call with `B > 1`: `min(B, W)`
`P`	Resident worker threads already owned by the current process

B controls numerical partitioning, ordered reduction, and failure placement. Reducing W does not renumber the B logical portions. Each of the J runners can process several portions and receives a stable scratch slot. The legacy one-thread, small, and empty paths execute directly. A prepared B > 1, W = 1 call made outside a worker still queues one runner. Nested prepared work executes all B portions locally with scratch slot zero and does not change queue counters.

The process pool keeps the largest worker capacity requested so far, so a later call can have P > W and P > J. This history does not increase that call's scratch allocation: planned scratch is sized from J, while partial results that preserve the numerical order remain sized from B or R. Separate concurrent calls own separate scratch and result buffers.

Native diagnostics expose cumulative process counters. worker_count is P; worker_start_events counts created resident workers; batches_submitted counts successfully committed queued calls; and tasks_submitted counts their queued runners, so a call with B > W adds J, not B. peak_queued_tasks is the largest observed queue depth. Direct calls add no batch or task count. These counters describe scheduling and are not a count of numerical rows, cells, or optimizer calls.

Memory budgets apply to one invocation. They check simultaneously live outputs, partials, prepared values, scratch slots, and binding-owned copies before work begins. They do not reserve process-wide capacity for other calls, and the resident pool's thread stacks are process resources rather than part of a kernel's memory_budget_bytes.

Parallelized workloads¶

The native implementation parallelizes independent blocks while preserving the sequential time recursions required by GAS and SCAR:

Workload	Parallel axis	Notes
Stochastic Student emission and gradient grids	observation rows	Reuses worker-local Student workspaces and the PPF cache
Equicorrelation Gaussian emission grids	observation rows	Reuses per-row normal-score sufficient statistics
Static multivariate row likelihoods and MLE	observation rows	Uses a prepared evaluator across optimizer calls
Gaussian and Student conditional sampling	generated rows	Reuses one conditional factorization when correlation is shared
Factor Student unconditional/conditional sampling	generated rows	Uses fixed Python draws and a compact native factor transform
Factor Gaussian likelihood and sampling	observation/generated rows	Reuses the immutable Woodbury operator; conditioning solves only `k*k`

GAS state updates and SCAR forward/backward filtering remain sequential over time. Increasing n_threads therefore accelerates the independent emission, row, or grid work, not the recurrence itself. Small workloads use a sequential fast path even when a larger value is requested, because thread scheduling would cost more than the kernel.

Multivariate sample, predict, and sample_conditional preserve the explicit n_threads setting through the model and top-level api entry points. Equicorr and stochastic Student batch methods preserve it for every block as well. Empty or fully fixed given mappings do not bypass thread validation. This setting controls observation sampling; it does not change the mathematical order of GAS or SCAR state updates.

Correctness and thread safety¶

Each submitted batch captures the caller's C floating-point environment and applies it on workers before numerical work. Worker defaults therefore do not change the arithmetic contract between serial and parallel execution.

Parallel kernels use a stable block partition and deterministic result placement. With identical inputs and random draws, the tested row, grid, and conditional-sampling paths agree between n_threads=1 and the parallel modes. If several rows fail, the parallel implementation reports the same smallest failure_index as the sequential implementation.

Random numbers are generated by the caller-facing Python layer and passed to native kernels as fixed draws. Use a fresh seeded generator when comparing thread counts:

one = student.sample_conditional(
    10_000,
    r=6.0,
    given={0: 0.4},
    rng=np.random.default_rng(17),
    n_threads=1,
)
many = student.sample_conditional(
    10_000,
    r=6.0,
    given={0: 0.4},
    rng=np.random.default_rng(17),
    n_threads=4,
)
np.testing.assert_array_equal(one, many)

Independent model instances can run concurrently. Mutating operations on the same model instance, such as simultaneous fits, are serialized by an instance-level lock. A prepared native evaluator also protects its mutable workspace; sharing one evaluator is safe but does not provide parallel objective calls. Use one model/evaluator per independent task when concurrent fits are desired.

The native pool belongs to the current process. Spawned or forked workers do not reuse a parent's pool. The pool also supports nested native calls through a worker-local sequential fallback, avoiding pool starvation and deadlock.

Embedded CPython subinterpreters are intentionally not supported. The native module declares multiple_interpreters::not_supported() and import from a subinterpreter must fail immediately. Ordinary Python threads and independent worker processes remain supported under the ownership rules above.

Independent fits and rolling windows¶

Use fit_independent when datasets, bootstrap replicas, model prototypes, or starting points are independent:

from pyscarcopula.contrib import fit_independent

batch = fit_independent(
    StochasticStudentCopula(d=d, R=R),
    bootstrap_samples,
    method="scar-tm-ou",
    fit_kwargs={"maxiter": 100},
    n_jobs=4,
)

print(batch.diagnostics)

In a Python script on a spawn-based platform such as Windows, place calls that create worker processes under the standard main guard:

if __name__ == "__main__":
    batch = fit_independent(model, datasets, n_jobs=4)

Each task reconstructs and owns its model, fit state, and prepared evaluator. With n_jobs > 1, omitting n_threads gives one native thread per process. Nested process/thread parallelism requires an explicit value:

batch = fit_independent(
    model,
    datasets,
    n_jobs=2,
    n_threads=4,
)

The same policy applies to risk_metrics. Per-window SeedSequence children preserve each window's random draws when process chunking changes. In an optimized portfolio, sequential execution carries the preceding window's weights into the next optimization. Each process chunk starts with equal weights and carries optimized weights only within that chunk. Changing n_jobs or chunk boundaries can therefore change optimized weights, VaR, and CVaR even with the same seed. Diagnostics report n_jobs, n_threads, the multiprocessing start method, whether nested parallelism was enabled, and the per-task ownership policy.

Avoid choosing n_jobs * n_threads substantially above the CPUs available to the job. Start with one parallelism level: native threads for one large fit, or processes with n_threads=1 for many independent fits.

Memory and dimensional limits¶

Parallel threads do not change the asymptotic representation of a model.

EquicorrGaussianCopula row and emission calculations use the scalar equicorrelation structure and avoid a dense correlation matrix in their hot path. prepare_sufficient_statistics accepts ndarray, memmap, and streamed blocks and stores only two O(T) vectors. Fixed dimension-tile reduction is deterministic across thread counts, and its default n_threads=1 does not initialize the native pool. pdf_and_grad_on_grid_batches bounds the (T,K) output, while memory_budget_bytes rejects an oversized monolithic output before native allocation. Unconditional sample_at_parameter_batches uses O(batch_rows*d) output memory and the structural equicorrelation eigenspaces for negative as well as positive correlation. Fitted sample_batches and predict_batches retain the sequential GAS update and SCAR path/posterior semantics while bounding each materialized output block. Monolithic sampling methods accept the same pre-allocation budget contract.
StochasticStudentCopula with corr_mode="fixed", "shrinkage", or "cholesky" still uses a dense correlation representation with O(d^2) storage/factorization costs.
Its corr_mode="factor" adapter stores O(d*k + k^2) state and routes static row and tiled latent-grid evaluation through the factor kernels. Supplied loadings or explicit two-stage initialization never build a dense covariance matrix. Static MLE, GAS, and SCAR-TM-OU likelihood consume the same immutable operator. SCAR emission selects independent cells or dimension tiles according to the workload. Unconditional and conditional Student generation uses the same operator; conditioning builds only a k*k factor system. Row batches bound the n*d output, and fixed seeds give identical results across thread counts. Static MLE additionally supports factor_estimation="joint" under the explicit factor_joint_max_params guard. Analytical loading gradients use fixed reduction blocks whose partition is independent of n_threads, so one- and multi-thread results are exact. The reduction workspace is a bounded constant multiple of d*k; no d*d gradient is formed. Joint loading estimation is available only for static MLE; GAS and SCAR use two-stage loadings.
GaussianCopula(corr_mode="factor") composes the same operator for native static likelihood, compact MLE, normal sampling, bounded batches, and exact conditioning. Tiled two-stage estimation, persistence, rolling workers, and the rank-dimensional Rosenblatt transform preserve the compact representation. A parallel Factor Rosenblatt call queues one row batch. Its shared-preparation path keeps four rank values per active scratch slot, plus alignment padding, and scratch is planned from J rather than the retained pool size P. The dense Gaussian mode remains the default.
The independent FactorCorrelation representation stores O(d*k + k^2) values and exposes prepared Woodbury matrix products, solves, quadratic forms, log determinants, and normal sampling. Its row kernels have an unconditional n_threads=1 default and accept explicit native row parallelism. Bounded sampling prevents an oversized monolithic output allocation.
FactorStudentEvaluator composes the factor operator with immutable observations for static Student row log densities and analytical derivatives with respect to df. Rows are the parallel axis; one active worker uses O(d) quantile/derivative/solve workspace. The default remains one thread and does not initialize the native pool.
Its tiled grid path does not retain an O(T*K*d) PPF table. It accumulates fixed dimension tiles directly into O(k) Woodbury and marginal summaries. Independent (row, df) cells are normally parallel; for a small grid and very large dimension, fixed dimension tiles are parallel and merged in deterministic order. evaluate_grid_batches bounds the (T,K) output and includes native/Python coexistence plus partial reductions in its peak-memory check.
A cached dynamic Student PPF table can require O(T * K * d) memory. It is capped at 256 MiB; when the values table is skipped, exact and asymptotic quantile paths preserve correctness with different performance tradeoffs.
Conditional sampling returns an O(n * d) result, so use application-level batches when the requested output itself is large.

Consequently, CPU threading alone does not make dense Student correlation modes suitable for d >> 10^4. Factor mode removes the dense correlation and PPF-cache limits for initialization, row/grid evaluation, MLE, GAS, and SCAR fitting. A large n_threads value by itself does not change the representation; corr_mode="factor" must be selected explicitly.

The compact Equicorr prepared object is consumed directly by static MLE, row/grid evaluation, GAS, and SCAR-TM-OU. These paths retain only the two statistic vectors after preparation. Streaming prediction/sampling output is a separate concern because the returned samples themselves require O(n * d) storage.

Native linear algebra¶

The extension uses a dependency-free C++17 linear-algebra layer for dense matrix-vector products, Cholesky factorization, SPD solves, and triangular products. It does not create Eigen, BLAS, OpenMP, or other external thread pools.

The runtime uses portable compiler-vectorizable reductions only for kernels whose end-to-end benchmarks pass the performance gate and retains the scalar backend for small kernels and regression comparisons. This keeps native threading under the explicit n_threads contract and avoids hidden BLAS oversubscription.

For optimizer and approximation controls, see Performance Tuning. For process-level helper signatures, see the Contrib API.