It achieves a better score 99 percent, then a fully connected neural network or SVM based on MNIST. Also, fast inference in us micro seconds.

A general N → M (f64 → f64) function approximator that reaches the accuracy of a well‑designed fully‑connected neural network on typical problems, but trains and infers much faster — a single convex least‑squares solve instead of iterative gradient descent.

It is not a neural network. There are no neurons, activations, layers, or backpropagation, and no neural‑network‑mimicking random features. Training has no local minima.

== Accuracy / timing on held-out test data ==
trig4       4->2   terms=131  train=0.011s  infer=0.52us/sample  R2=1.00000  RMSE=7.5e-06
friedman1   10->1  terms=776  train=0.413s  infer=1.54us/sample  R2=1.00000  RMSE=0.0107
multi6      6->3   terms=286  train=0.046s  infer=0.76us/sample  R2=1.00000  RMSE=6.6e-06

The fit is one convex least-squares problem, but it is solved two ways depending on the number of retained terms F (auto-selected via dense_max):

• Dense path (F ≤ dense_max, default 1200): the normal equations ΦᵀΦ are assembled in parallel across all CPU cores and solved with a Cholesky factorization shared across the M outputs. Exact and extremely fast for the typical low/moderate-dimensional case.
• Scalable path (F > dense_max): a matrix-free, Jacobi-preconditioned block Conjugate Gradient solver. It never forms an F×F matrix (so memory is O(F)), applies ΦᵀΦ matrix-free and parallelized across all cores, and reaches the same unique convex minimizer. The design matrix is cached in f64 when it fits a memory budget (for a low residual floor; coefficients are accumulated in f64 as well).

Relevant Config knobs: dense_max, threads (0 = auto), cg_iter, cg_tol.

== Accuracy / timing on held-out test data ==
trig4       4->2   terms=131    train=0.004s  infer=0.31us  R2=1.00000
friedman1   10->1  terms=776    train=0.096s  infer=1.52us  R2=1.00000
multi6      6->3   terms=286    train=0.011s  infer=0.78us  R2=1.00000

== Dataset-size scaling (friedman1, 10->1) ==
rows=10000   train=0.10s
rows=50000   train=0.31s
rows=200000  train=1.80s

== High dimensional (1000 inputs, 8 relevant, 992 noise -> 1) ==
rows=20000   terms=10001  train=1.38s  test-R2=1.00000
rows=50000   terms=10001  train=1.87s  test-R2=1.00000

At 1000 input dimensions the model fits to R² = 1.0 in ~1.4 s and correctly ignores the 992 irrelevant inputs (the ridge penalty drives their coefficients to ~0; the R² is on held-out data, so this is genuine generalization, not overfitting). It scales further (2000-D fits to R² ≈ 0.998).

The number of terms grows combinatorially with the interaction order, so at very high dimensions full pairwise (order 2) is impractical (≈7.5M terms at 1000-D). For high-dimensional inputs use the additive regime (max_order = 1), which captures arbitrary per-dimension nonlinearity and is what the 1000-D benchmark above uses. Lower-dimensional problems use order 2–3 for interactions.

The same approximator is applied to the real MNIST handwritten-digit dataset (60000 train / 10000 test, 28×28 grayscale, 10 classes). Classification is cast as multiclass least-squares classification: the 10 classes become one-hot targets in R^10, the convex regressor fits all 10 outputs in a single solve, and a test image is labeled by argmax over its 10 (de-standardized) predicted outputs. Still no neural network, no gradient descent, no random features.

Because raw pixels are highly correlated and a full pairwise expansion over 784 inputs is infeasible, each image goes through two fixed, deterministic preprocessing steps (no learning, no randomness — so this stays non-neural and free of random features):

1. Moment-based deskewing — estimate each digit's slant from its second-order pixel moments and undo the shear (bilinear resampling), removing most of the handwriting-slant variation.
2. 2D DCT — keep the low-frequency Kf×Kf coefficients. The DCT is the classic image decorrelator (the basis of JPEG); the coefficients are near-uncorrelated and few, which conditions the convex solve and makes an all-pairs degree-2 (polynomial-kernel-like) expansion feasible.

Measured on a 24-core CPU (./fga mnist):

deskew (moment-based shear)                (0.44s)
2D-DCT 10x10 -> 100 decorrelated features  (0.14s)
DCT(10x10) additive     terms=601    train=0.21s    infer=2.7us    TEST-ACC=0.9589  (err 4.11%)
DCT(10x10) + pairwise   terms=5551   train=30.3s    infer=26us     TEST-ACC=0.9833  (err 1.67%)
DCT(10x10) + graded234  terms=24667  train=98.2s    infer=118us    TEST-ACC=0.9885  (err 1.15%)
2D-DCT 12x12 -> 144 decorrelated features  (0.17s)
DCT(12x12) additive     terms=865    train=0.61s    infer=4.0us    TEST-ACC=0.9608  (err 3.92%)
DCT(12x12) + pairwise   terms=11161  train=82.0s    infer=51us     TEST-ACC=0.9854  (err 1.46%)
DCT(12x12) + graded234  terms=30277  train=124.4s   infer=144us    TEST-ACC=0.9886  (err 1.14%)
2D-DCT 13x13 -> 169 decorrelated features  (0.19s)
DCT(13x13) additive     terms=1015   train=1.36s    infer=4.6us    TEST-ACC=0.9618  (err 3.82%)
DCT(13x13) + pairwise   terms=15211  train=135.4s   infer=70us     TEST-ACC=0.9858  (err 1.42%)
DCT(13x13) + graded234  terms=34327  train=145.1s   infer=163us    TEST-ACC=0.9896  (err 1.04%)
2D-DCT 14x14 -> 196 decorrelated features  (0.20s)
DCT(14x14) additive     terms=1177   train=4.14s    infer=5.4us    TEST-ACC=0.9620  (err 3.80%)
DCT(14x14) + pairwise   terms=20287  train=201.0s   infer=91us     TEST-ACC=0.9854  (err 1.46%)

Adding the degree-2 interactions lifts accuracy from 95.89% → 98.33% (10×10). Larger DCT blocks keep helping with diminishing returns (13×13 pairwise reaches 98.58%). Adding graded higher-order interactions — on top of the all-pairs degree-2 base, the top-48 features (by a cheap relevance screen) get degree-3 (triple) cross terms and the top-16 features get degree-4 (quadruple) cross terms — lifts every block further, with 13×13 + graded234 reaching the best 98.96% in ~145 s, and even 10×10 + graded234 hitting 98.85% in ~98 s. The additive model alone trains in 0.2 s. (Without deskewing the two degree-2 models score 93.3% / 97.4%, so deskewing alone is worth ~+1 point.)

Reading: this method matches and slightly exceeds the best single fully- connected neural network (98.33–98.96% vs ~98.4%) and beats classical baselines like k-NN and a linear classifier — while training as one convex least-squares solve on a CPU in tens of seconds (no GPU, no epochs, no backprop, no local minima) and inferring in ~30 µs/sample. With screened graded degree-2/3/4 interactions it reaches 98.96%, matching a tuned polynomial/RBF kernel SVM, while convolutional-network ensembles remain ahead in absolute accuracy because they exploit image structure this general-purpose approximator does not. The point is not to beat a CNN, but to show a genuinely different, convex, non-neural approximator reaching (single-)neural-network-class accuracy on a real task, fast.

The honest caveats: (1) deskewing and DCT are image-specific preprocessing (the generic approximator itself is unchanged); (2) least-squares one-hot classification is not softmax — it is the classic linear-discriminant view; (3) we report a single test evaluation with hand-set hyperparameters (Kf, degrees, ridge), tuned only on the training split. Numbers for the comparison rows are well-known literature values.

mkdir -p data && cd data
for f in train-images-idx3-ubyte train-labels-idx1-ubyte \
         t10k-images-idx3-ubyte  t10k-labels-idx1-ubyte; do
  curl -sSL https://storage.googleapis.com/cvdf-datasets/mnist/$f.gz | gunzip > $f
done

Then odin build src -out:fga -o:speed && ./fga mnist. (Set DO_RAW_ADDITIVE to true in src/mnist.odin to also run the slow ~2-min raw-pixel additive baseline, which tops out at ~89%.)

The target is approximated by a truncated functional‑ANOVA (a.k.a. sparse polynomial‑chaos) expansion in a tensor‑product Chebyshev basis:

f(x) ≈ c0                                   (constant)
     + Σ_d   Σ_p      c[d,p]·T_p(x_d)               (univariate / additive part)
     + Σ_{d<e} Σ_{p,q} c[d,e,p,q]·T_p(x_d)·T_q(x_e)  (pairwise interactions)
     + (optional triple and quadruple interactions)

where T_k is the degree‑k Chebyshev polynomial of the first kind (extremely well conditioned on [-1,1]). Each input is rescaled to [-1,1]; outputs are standardized for conditioning.

• Linear in the coefficients ⇒ fitting is one globally convex ridge least‑squares problem. There are no local minima, so the fit is reliable and reproducible — unlike alternating tensor (TT/ALS) or neural methods.
• One solve, no iterations. We assemble the small dense normal equations A = ΦᵀΦ, rhs = ΦᵀY in a single streaming pass over the data, normalize the feature columns to unit scale (so a single small ridge regularizes uniformly), and solve with a Cholesky factorization shared across all M outputs.
• Curse of dimensionality is tamed by truncation. Only low‑order interactions (default: order ≤ 2) with bounded total degree are kept. This natively represents additive structure and low‑order interactions — the structure most real functions actually have — while keeping the number of terms small.
• Inference is microseconds: evaluate the per‑dimension Chebyshev tables once, then a sum of sparse products.

• Long‑range additive structure such as y = x₀ + x₅ (the worst case for low‑rank tensor‑train ALS, which gets trapped in local minima) is fit to R² = 1.0 here, because additive terms are first‑class basis functions.
• Interaction order is a principled, honest tuning knob. A genuine 3‑way product x₀·x₁·x₂ is approximated to R² ≈ 0.90 at order 2 and R² = 1.0 at order 3.

Polynomial‑chaos expansions and the functional‑ANOVA decomposition from uncertainty quantification and surrogate modelling (Wiener; Sobol’; Blatman & Sudret, sparse PCE). Distinct from neural networks and from tensor‑train ALS.

import fga "src"

cfg := fga.default_config()        // deg1=10, deg2=6, deg3=3, max_order=2, ridge=1e-6
m   := fga.train(X, Y, rows, n_in, n_out, cfg)
defer fga.destroy_model(&m)

s   := fga.make_predict_scratch(&m)
defer fga.destroy_predict_scratch(&s)
out := make([]f64, n_out); defer delete(out)
fga.predict(&m, x, out, &s)        // x: []f64 length n_in -> out length n_out

X is rows*n_in row‑major, Y is rows*n_out row‑major.

For high-dimensional inputs set max_order = 1 (additive). Increase max_order/degrees for stronger high‑order interactions (more terms, slower); decrease them for very high‑dimensional inputs. For max_order ≥ 3, set triple_dims / quad_dims to screened subsets of informative features to keep the C(n,3) / C(n,4) term counts tractable (the MNIST demo screens triples to the top 48 and quadruples to the top 16 features).

odin build src -out:fga -o:speed
./fga

Requires the Odin compiler. The demo (src/main.odin) reports accuracy and timing on three synthetic benchmarks plus a training‑time scaling test.

• src/approx.odin — model, Chebyshev basis, sparse‑term evaluation, predict, metrics.
• src/train.odin — term‑set generation and the convex fit (dense or CG path).
• src/solver.odin — parallel fork‑join, parallel dense normal‑equations assembly, and the matrix‑free parallel Conjugate‑Gradient solver.
• src/linalg.odin — Cholesky factor/solve and the multi‑RHS SPD solver.
• src/main.odin — benchmark / demo harness (incl. the 1000‑D test).
• src/mnist.odin — real-MNIST benchmark (2D-DCT features + one-hot least-squares classification); run with ./fga mnist.

• Training: one pass to assemble the normal equations, O(rows · F²) where F is the (small) number of retained terms, plus one O(F³) Cholesky. Linear in the number of samples.
• Inference: O(F) per sample after O(n_in · deg_max) basis evaluation.

The Discrete Cosine Transform (DCT) rewrites a signal as a sum of cosine waves of increasing frequency. It is a change of basis: the same information, expressed in terms of "how much of each frequency is present" instead of raw sample values. It is the transform at the heart of JPEG and MP3.

For a length-N signal x[0..N-1], the DCT-II produces N coefficients:

X[k] = Σ_{n=0}^{N-1}  x[n] · cos( π/N · (n + ½) · k ) ,     k = 0 … N-1

• X[0] (k=0) is just the sum/average — the DC / lowest-frequency component (the cosine of frequency 0 is the constant 1).
• Increasing k measures faster and faster oscillations — higher frequencies.
• The basis vectors cos(π/N·(n+½)·k) are mutually orthogonal, so this is an exact, invertible, energy-preserving change of basis (like rotating the coordinate axes).

Intuition. Each coefficient answers "how much does the signal look like this cosine wave?" A smooth signal is mostly low-frequency, so its energy piles up in the first few coefficients and the high-frequency ones are ≈0. This is called energy compaction, and it is why you can throw away most high-frequency coefficients with little loss (JPEG does exactly this).

1. Energy compaction → dimensionality reduction. Natural images are smooth, so almost all of a digit's "shape energy" lives in the low-frequency coefficients. Keeping only the top-left Kf×Kf block (e.g. 10×10 = 100 of the 784) discards little useful information. → 784 inputs become ~100.

2. Decorrelation → better conditioning. The DCT closely approximates the Karhunen–Loève transform (a.k.a. PCA) for the kind of correlations natural images have. The output coefficients are therefore nearly uncorrelated, even though the input pixels are strongly correlated. Uncorrelated features make the regression's normal matrix ΦᵀΦ nearly diagonal — well-conditioned — so the convex solver converges fast and stably.

An image is 2-D, so we apply the 1-D DCT twice: once along every row, then once along every column of the result (the order doesn't matter). Mathematically:

F[u,v] = Σ_r Σ_c  image[r,c] · cos(π/N·(r+½)·u) · cos(π/N·(c+½)·v)

Because it separates into two 1-D passes, the cost is just 2·N cosines per output coefficient instead of N² — cheap. In the code (dct_features) we precompute the Kf×28 matrix of cosine values once, then for each image do a rows-pass into a small Kf×28 buffer and a columns-pass into the Kf×Kf output. F[0,0] is the overall brightness; F[u,v] for small u,v captures the coarse stroke layout; large u,v would capture fine detail/noise, which we drop.

DCT = "describe the image by its coarse-to-fine cosine ‘ingredients’." Keep the coarse ingredients (low frequencies): you get a small, decorrelated, information- rich feature vector that a convex model can fit easily.

The headline result (98.33–98.96% test accuracy, trained as one convex solve on a CPU) comes from stacking several deterministic, principled ideas. None of them is a neural network, gradient descent, or a random feature. In order of impact:

1. Moment-based deskewing (deskew_images). The biggest nuisance variation in MNIST is slant. We measure each digit's slant directly from its image moments: the centre of mass (m_r, m_c), the vertical spread var_r, and the row/column covariance cov. The slant is α = cov / var_r. We then apply the inverse horizontal shear (column ← column − α·row), recentre on the centre of mass, and resample with bilinear interpolation. This makes a "7" written upright and a slanted "7" look the same to the model — variation a convex classifier cannot otherwise remove. Worth ≈ +1 point (97.4 → 98.3 at 10×10).

2. 2D-DCT feature reduction (Appendix A). Converts 784 correlated pixels into ~100 decorrelated, low-frequency coefficients. Two wins at once: (a) it conditions the convex solve (uncorrelated features → near-diagonal ΦᵀΦ → fast CG); (b) it shrinks dimensionality enough that an all-pairs degree-2 expansion is feasible (C(100,2)=4 950 pairs instead of C(784,2)=306 936).

3. One-hot least-squares classification + de-standardized argmax. The 10 classes become one-hot target vectors in R¹⁰; the regressor fits all 10 outputs in a single convex solve; a test image is labeled by argmax over the 10 predicted outputs (after un-standardizing, so class priors enter correctly). This is the classic linear-discriminant / least-squares classifier view — convex, no iteration beyond the one solve.

4. Degree-2 (pairwise) interactions = a polynomial classifier. Univariate terms alone (a generalized additive model) reach 95.89%. Adding all-pairs products of the DCT features makes the model a degree-2 polynomial in feature space — the explicit primal form of a polynomial-kernel SVM — and lifts accuracy to 98.33% (10×10). This is the single largest modelling gain.

4b. Graded higher-order interactions (degree 2/3/4), screened by feature rank. A full degree-3 polynomial would add C(fdim,3) cubic cross terms (≈162k at 10×10, ≈794k at 13×13) and degree-4 even more — too many to fit fast. Instead each feature is ranked by a cheap, label-driven relevance score (the sum over classes of its squared correlation with the one-hot target), and higher orders are reserved for the most informative features: every feature keeps its degree-2 (all-pairs) terms, the top 48 features additionally get degree-3 (triple) terms (C(48,3) = 17 296), and the top 16 features additionally get degree-4 (quadruple) terms (C(16,4) = 1 820). This keeps the model small and well-conditioned, yet lifts accuracy another notch: 98.85% (10×10) up to the best 98.96% (13×13), all still in under ~150 s of CPU training.

5. Chebyshev basis for per-feature nonlinearity. Each feature's nonlinearity is expanded in Chebyshev polynomials T_0…T_6, which are near-optimally conditioned on [-1,1] (unlike raw powers x, x², …, which are nearly collinear and numerically unstable). This lets us use high per-feature degree without blowing up the condition number.

6. Per-feature (column-normalized) ridge. Instead of a flat λ·I, the ridge is scaled by each feature's own energy (R = diag(λ · diag(ΦᵀΦ))), so regularization acts uniformly across features of very different magnitudes (DC coefficient vs. a high-frequency one). Keeps the solve well-posed and lets weak features shrink cleanly.

7. A deliberately larger ridge for the iterative path. Bumping λ to 3e-3 improves the condition number, so CG needs far fewer iterations — and it doubles as healthy regularization for the classifier.

8. f64 design-matrix cache + early iteration cap. The CG mat-vec is memory-bandwidth-bound, so the design matrix is cached once (rather than recomputed every iteration) in 64-bit floats, giving a very low residual floor. Classification accuracy saturates early regardless — ~80 iterations give the same accuracy as 300 — so we cap iterations low and finish several times sooner.

9. Block Conjugate Gradient over all 10 outputs. All 10 one-hot columns are solved simultaneously, sharing the one expensive ΦᵀΦ·p mat-vec per iteration, so 10-output training costs essentially the same as 1-output.

10. Multithreaded, matrix-free solve. The mat-vec (and the dense assembly) are parallelized across all CPU cores with a lock-free fork-join (per-thread partials, then a reduce), and CG never forms the F×F matrix — keeping memory at O(F) and the whole fit at tens of seconds on a CPU.

Everything above is fixed and deterministic: same data ⇒ same model, every time.

How long does it take to reach this accuracy on a CPU (no GPU)? The figures below are representative single-machine, multicore-CPU numbers; absolute times vary with hardware, BLAS, and implementation, but the orders of magnitude are the point. This method's numbers are measured on a 24-core CPU; the others are typical values for standard CPU implementations (e.g. scikit-learn / libsvm) reaching the listed accuracy.

Takeaways.

• Accuracy class. At 98.33–98.96% this method matches/edges out a well-tuned single fully-connected neural network (~98.4%), beats k-NN, logistic regression and random forests, and trails only kernel SVMs and convolutional nets.
• Training speed. It reaches that accuracy in one convex least-squares solve — seconds to ~2 minutes on a CPU — versus minutes to tens of minutes for an MLP (which must run many SGD epochs) or a kernel SVM (whose training is quadratic in the 60 000 samples). The additive model (95.9%) trains in under a second.
• No GPU, no epochs, no tuning loop. There is no learning-rate schedule, no early-stopping, no random initialization, and no local minima — the same data always yields the same model. An MLP of comparable accuracy needs epoch sweeps and hyperparameter/seed search to match it.
• Inference. Predicting is a short sparse dot product: 3–60 µs/sample, comparable to an MLP and far faster than k-NN (which rescans the training set) or a kernel SVM (which sums over many support vectors).

The honest counterpoint: convolutional networks still win on raw accuracy (99%+), because convolutions encode image structure (locality, translation invariance) that this general-purpose approximator does not. The achievement here is reaching single-MLP-class accuracy with a completely different, convex, non-neural method that trains by solving one linear system.

MIT Open Source License

Best regards,
Joao Carvalho