A fast, header-only C library for Singular Spectrum Analysis with Python bindings, FFT-accelerated operations, and Intel MKL optimization.

17-30× faster than Rssa (the gold-standard R implementation).

Singular Spectrum Analysis decomposes time series into:

• Trends - Long-term movements
• Periodic components - Cycles, seasonality
• Noise - Random fluctuations

Unlike Fourier analysis, SSA is non-parametric and adapts to your data's structure. Ideal for financial analysis, signal processing, and anywhere you need clean signal extraction.

SSA recovers clean signal from heavy noise (SNR: 4.1 → 26.1 dB).

Ensemble SSA extracts smooth trend from volatile Bitcoin prices (r = 0.987).

2-year ahead forecast capturing trend and 12-month seasonality (RMSE: 47.5).

Benchmarked with benchmark_vs_rssa.py on Intel Core i9-14900KF, MKL 2025.0.3:

• Components: k = 50
• Window length: L = N/4 (standard SSA choice)
• Signal: Synthetic (trend + periodic + noise)

Summary: Average 24× speedup, peak 30.5× at N=10,000. Throughput: 625 SSA/sec at N=1000.

Rssa uses PROPACK (Lanczos bidiagonalization). SSA-Opt uses randomized SVD with MKL.

Using mkl_config.h for hybrid CPU optimization (P-core affinity, thread tuning):

The mkl_config.h header auto-detects hybrid CPUs (Intel 12th-14th gen) and pins threads to P-cores only, avoiding E-core slowdown.

• 3 Decomposition Methods - Sequential power iteration, block power iteration, and randomized SVD
• FFT-Accelerated - O(N log N) Hankel matrix-vector products via convolution theorem
• R2C Optimization - Real-to-complex FFT exploits Hermitian symmetry (50% memory savings)
• Malloc-Free Hot Path - prepare() + decompose_randomized() for zero-allocation streaming

• Grouped Reconstruction - Combine arbitrary component subsets
• W-Correlation Matrix - Fast DSYRK-based weighted correlation for component grouping
• Automatic Pair Detection - Find sine/cosine pairs via singular value + W-correlation matching
• Variance Explained - Per-component and cumulative variance statistics
• Component Statistics - Scree plot data, gap detection, automatic rank selection

• R-Forecast (LRF) - Linear Recurrent Formula for extrapolation
• V-Forecast - Alternative vector forecast method
• Verticality Check - Stability indicator for forecast reliability

• ESPRIT - Frequency/period estimation from eigenvectors
• Cadzow Iterations - Iterative rank reduction for enhanced denoising
• Gap Filling - Missing value imputation (iterative + simple methods)
• MSSA - Multivariate SSA for multi-channel/correlated series

• MKL-Optimized - Intel MKL for FFT, BLAS, LAPACK, and RNG
• Batched FFT - Process multiple vectors in single MKL call
• Cached FFTs - Pre-computed eigenvector FFTs for fast reconstruction
• OpenMP Ready - Parallel-friendly design

• Header-Only - Single #include "ssa_opt.h"
• Python Bindings - Full ctypes wrapper with NumPy integration
• Minimal Dependencies - Just Intel MKL (or compatible BLAS/LAPACK)

from ssa_wrapper import SSA # Load your data prices = np.array([...]) # Your time series # Decompose ssa = SSA(prices, L=250) ssa.decompose(k=30) # Uses fast randomized SVD # Extract components trend = ssa.reconstruct([0]) cycle = ssa.reconstruct([1, 2]) # Forecast forecast = ssa.forecast([0, 1, 2], n_forecast=50) # Analysis print(f"Variance explained: {ssa.variance_explained(0, 2):.1%}") pairs = ssa.find_periodic_pairs()

#define SSA_OPT_IMPLEMENTATION #include "ssa_opt_r2c.h" #include "mkl_config.h" int main() { // Initialize MKL (call once at startup) mkl_config_ssa_full(1); // 1 = verbose // Setup SSA_Opt ssa = {0}; ssa_opt_init(&ssa, signal, N, L); // Decompose (randomized is fastest) ssa_opt_decompose_randomized(&ssa, k, 8); // Reconstruct trend int trend_group[] = {0}; double* trend = malloc(N * sizeof(double)); ssa_opt_reconstruct(&ssa, trend_group, 1, trend); // Cleanup ssa_opt_free(&ssa); free(trend); }

Access results directly via the struct:

ssa->sigma[i] // Singular values (descending order) ssa->eigenvalues[i] // Squared singular values (variance explained) ssa->U[i * L] // Left singular vector i (length L) ssa->V[i * K] // Right singular vector i (length K)

SSA embeds a 1D time series into a 2D trajectory matrix. Given signal x = [x₀, x₁, ..., x_{N-1}] and window length L, we construct an L × K Hankel matrix where K = N - L + 1:

 Column 0 Column 1 Column 2 ... Column K-1 ┌─────────────────────────────────────────────────┐ Row 0 │ x₀ x₁ x₂ ... x_{K-1} │ Row 1 │ x₁ x₂ x₃ ... x_K │ Row 2 │ x₂ x₃ x₄ ... x_{K+1} │ ⋮ │ ⋮ ⋮ ⋮ ⋱ ⋮ │ Row L-1│ x_{L-1} x_L x_{L+1} ... x_{N-1} │ └─────────────────────────────────────────────────┘ H[i,j] = x[i+j] (constant along anti-diagonals) 

This matrix is never explicitly formed - it would require O(L × K) ≈ O(N²) memory.

We compute the top-k singular triplets of H:

H ≈ σ₀·u₀·v₀ᵀ + σ₁·u₁·v₁ᵀ + ... + σₖ₋₁·uₖ₋₁·vₖ₋₁ᵀ 

Each triplet (σᵢ, uᵢ, vᵢ) represents one component. Larger σ = more signal energy.

Each component forms a rank-1 matrix, then averaged back to 1D along anti-diagonals:

t=0: x̃[0] = X[0,0] t=1: x̃[1] = (X[0,1] + X[1,0]) / 2 t=2: x̃[2] = (X[0,2] + X[1,1] + X[2,0]) / 3 ... 

SSA signals satisfy a linear recurrence. The LRF extracts coefficients from left singular vectors to predict future values:

x[n] = a₁·x[n-1] + a₂·x[n-2] + ... + a_{L-1}·x[n-L+1] 

Hankel matrix-vector multiplication equals convolution:

Forward: y = H·v Adjoint: z = Hᵀ·u y[i] = Σⱼ x[i+j]·v[j] z[j] = Σᵢ x[i+j]·u[i] = (x ∗ v)[i] = (x ∗ flip(u))[L-1+j] ┌─────────────────────────────────────────────────────────┐ │ conv(a,b) = IFFT( FFT(a) ⊙ FFT(b) ) │ │ │ │ Complexity: O(N²) direct → O(N log N) with FFT │ └─────────────────────────────────────────────────────────┘ Pre-compute FFT(signal) once at init. Each matvec: 1 FFT + 1 pointwise multiply + 1 IFFT 

Process b vectors simultaneously instead of sequentially:

Sequential (b=1): Block (b=32): ┌───┐ ┌───────────────────┐ │ v │ → H·v → u → Hᵀ·u → v │ v₀ v₁ v₂ ... v₃₁ │ = V_block └───┘ └───────────────────┘ ↓ × k components H·V_block (32 FFTs batched) × 100 iterations ↓ = 6400 individual FFTs QR(U_block) ↓ Hᵀ·U_block (32 FFTs batched) ↓ = 200 batched calls total Benefits: • MKL batched FFT: single function call for 32 transforms • Memory locality: V_block columns contiguous in memory • BLAS3: QR uses GEMM instead of GEMV 

Every iteration: vs. Every 5 iterations: iter 1: QR iter 1: (skip) iter 2: QR iter 2: (skip) iter 3: QR iter 3: (skip) iter 4: QR iter 4: (skip) iter 5: QR iter 5: QR ← ⋮ ⋮ 100 QR calls 20 QR calls (5× reduction) Stability maintained by: • Final QR before Rayleigh-Ritz extraction • Deflation orthogonalization against previous blocks each iteration 

After block iteration converges, vectors span the correct subspace but aren't individually optimal:

┌─────────────────────────────────────────────────────┐ │ 1. Project: M = U_blockᵀ · (H · V_block) [b × b] │ │ │ │ 2. Small SVD: M = Uₛ · Σ · Vₛᵀ [b × b] │ │ │ │ 3. Rotate: U_final = U_block · Uₛ │ │ V_final = V_block · Vₛ │ │ σ_final = diag(Σ) │ └─────────────────────────────────────────────────────┘ Cost: one 32×32 SVD (microseconds) → optimal singular values 

All multi-vector operations use BLAS Level 3:

GEMM achieves near-peak FLOPS due to cache reuse and SIMD.

For very large problems, randomized SVD is fastest:

1. Random projection: Y = H · Ω, Ω is K × (k+p) random 2. Orthogonalize: Q = qr(Y) 3. Project: B = Qᵀ · H 4. SVD of small matrix: B = U_B · Σ · Vᵀ 5. Recover: U = Q · U_B 

Complexity: O((k+p) × N log N) vs O(k × iterations × N log N) for power iteration.

Real-to-Complex FFT exploits conjugate symmetry:

• Output: N/2+1 complex values (vs N for C2C)
• 50% less memory for FFT buffers
• ~17% faster than complex-to-complex

W-correlation matrix requires FFTs of all component reconstructions. With caching:

Call ssa_opt_cache_ffts() once after decomposition if computing W-correlation multiple times.

# Build the shared library cd MKL mkdir build && cd build cmake .. -DUSE_MKL=ON cmake --build . --config Release # Copy ssa.dll/libssa.so to py/ folder # Then: cd ../py python -c "from ssa_wrapper import SSA; print('OK')"

#define SSA_OPT_IMPLEMENTATION #include "ssa_opt_r2c.h"

Link with MKL:

# Windows cl /O2 your_code.c /I"%MKLROOT%\include" mkl_rt.lib # Linux gcc -O3 your_code.c -I${MKLROOT}/include -lmkl_rt -lm

Run the benchmark suite:

cd py # Speed comparison vs Rssa python benchmark_vs_rssa.py --speed # Compare MKL configurations python benchmark_vs_rssa.py --mkl-config # Internal method comparison (seq/block/randomized) python benchmark_vs_rssa.py --internal # All benchmarks python benchmark_vs_rssa.py --all

Parameter sensitivity tests:

python ssa_parameter_tests.py --snr # Noise robustness python ssa_parameter_tests.py --window # L sensitivity python ssa_parameter_tests.py --components # k sweep python ssa_parameter_tests.py --all

• Memory: O(k × N) for storing singular vectors
• Not for streaming: Batch processing only
• Minimum size: N > 100 recommended

1. Golyandina, N., & Zhigljavsky, A. (2013). Singular Spectrum Analysis for Time Series. Springer.
2. Korobeynikov, A. (2010). Computation- and space-efficient implementation of SSA. Statistics and Its Interface, 3(3).
3. Halko, N., Martinsson, P. G., & Tropp, J. A. (2011). Finding structure with randomness. SIAM Review, 53(2).

GPL 3.0 - See LICENSE file.