Edgar Chraghyan · ISO Scientific · Advanced technical tutorial
Abstract
This tutorial develops contemporary machine-learning training from explicit statistical models, matrix calculus, numerical optimization and reproducible low-level implementations. The objective is to connect the mathematics to working C and C++17 programs rather than treating a training framework as a black box. It covers likelihood-based learning, gradient descent and Adam, numerical gradient checking, nonlinear classifiers, attention, self-supervised objectives, parameter-efficient adaptation and rigorous evaluation. The algorithms discussed here are standard, publicly described methods; the demonstrations are educational.
1. Statistical learning as an optimization problem
Let the examples be $D=\{(x_i,y_i)\}_{i=1}^{n}$, with $x_i\in\mathbb R^d$ and labels or targets $y_i$. A parameterized predictor $f_{\theta}:\mathbb R^d\rightarrow\mathcal Y$ is selected by minimizing an empirical objective:
$$ \widehat\theta=\arg\min_{\theta\in\mathbb R^p}\left\{\frac1n\sum_{i=1}^{n}\ell(f_\theta(x_i),y_i)+\lambda\Omega(\theta)\right\}. $$
The first term measures fit to observed data; $\Omega$ penalizes undesirable complexity. With an independent test distribution $P$, the quantity of interest is instead the population risk $R(\theta)=\mathbb E_{(x,y)\sim P}[\ell(f_\theta(x),y)]$. Low training loss is not sufficient evidence of low population risk. Random data splits must be constructed before tuning; time-series or grouped data need temporal/group-aware splitting instead of naïve random shuffling.
1.1 Linear regression and a complete derivation
For $X\in\mathbb R^{n\times d}$, $y\in\mathbb R^n$, and $\widehat y=Xw+b\mathbf1$, let
$$ J(w,b)=\frac1n\lVert Xw+b\mathbf1-y\rVert_2^2+\lambda\lVert w\rVert_2^2. $$
Write $r=Xw+b\mathbf1-y$. The differential is $d(r^Tr)=2r^Tdr$; since $dr=Xdw+\mathbf1db$, we obtain
$$ \nabla_wJ=\frac2nX^Tr+2\lambda w,\qquad \frac{\partial J}{\partial b}=\frac2n\mathbf1^Tr. $$
The intercept is not regularized in this definition. When $\lambda=0$ and the columns of $X$ are full rank (with an intercept handled explicitly), the normal equations lead to $w=(X^TX)^{-1}X^Ty$. Numerical software should solve the linear system by QR or SVD rather than explicitly forming this inverse. The Hessian in $w$ is $(2/n)X^TX+2\lambda I$, which is positive semidefinite and becomes positive definite when $\lambda>0$.
1.2 Likelihood and cross-entropy
For binary $y\in\{0,1\}$ and $p_\theta(x)=\sigma(z)$ with $z=w^Tx+b$ and $\sigma(z)=1/(1+e^{-z})$, the negative log-likelihood is
$$ \ell(z,y)=-y\log\sigma(z)-(1-y)\log(1-\sigma(z)). $$
Differentiate carefully: $d\sigma/dz=\sigma(1-\sigma)$ and thus $\partial\ell/\partial z=\sigma(z)-y$. It follows immediately that
$$ \nabla_wJ=\frac1n\sum_{i=1}^{n}(\sigma(z_i)-y_i)x_i,\qquad \partial_bJ=\frac1n\sum_{i=1}^{n}(\sigma(z_i)-y_i). $$
For stable evaluation, use $\ell(z,y)=\max(z,0)-zy+\log(1+e^{-|z|})$ instead of directly taking logarithms of probabilities close to zero.
2. Backpropagation from the chain rule
Consider a two-layer network with $W_1\in\mathbb R^{h\times d}$, $b_1\in\mathbb R^h$, $W_2\in\mathbb R^{k\times h}$, $b_2\in\mathbb R^k$:
$$ z_1=W_1x+b_1,\quad a_1=\phi(z_1),\quad z_2=W_2a_1+b_2,\quad L=\ell(z_2,y). $$
Define $\delta_2=\nabla_{z_2}L$. Reverse-mode differentiation gives
$$ \frac{\partial L}{\partial W_2}=\delta_2a_1^T,\quad \frac{\partial L}{\partial b_2}=\delta_2,\quad \delta_1=(W_2^T\delta_2)\odot\phi'(z_1),\quad \frac{\partial L}{\partial W_1}=\delta_1x^T. $$
The outer-product shape is important: $\delta_2a_1^T$ has dimensions $k\times h$. Broadcasting and transposition errors can silently corrupt training even when an implementation compiles. For a minibatch, average the per-example gradients (or use equivalent batched matrix operations).
2.1 Finite-difference gradient checks
Before using a complicated optimizer, compare each analytic gradient component with a centered finite difference:
$$ g_j^{FD}=\frac{L(\theta+\varepsilon e_j)-L(\theta-\varepsilon e_j)}{2\varepsilon},\quad \rho_j=\frac{|g_j-g_j^{FD}|}{\max(1,|g_j|,|g_j^{FD}|)}. $$
Start with double precision, $\varepsilon\approx10^{-5}$, small networks, and deterministic code. Too small an $\varepsilon$ amplifies floating-point cancellation; too large an $\varepsilon$ amplifies truncation error. Gradient checking validates an implementation locally; it does not establish generalization.
3. Optimizers and their limits
Full-batch gradient descent uses $\theta_{t+1}=\theta_t-\eta\nabla J(\theta_t)$. Stochastic gradient descent uses an unbiased or approximately unbiased minibatch estimate. With momentum,
$$ v_t=\mu v_{t-1}+g_t,\qquad \theta_t=\theta_{t-1}-\eta v_t. $$
Adam maintains first and second moments of $g_t$ and corrects initialization bias:
$$ m_t=\beta_1m_{t-1}+(1-\beta_1)g_t,\quad v_t=\beta_2v_{t-1}+(1-\beta_2)g_t^2, $$
$$ \widehat m_t=m_t/(1-\beta_1^t),\quad \widehat v_t=v_t/(1-\beta_2^t),\quad \theta_{t+1}=\theta_t-\eta\widehat m_t/(\sqrt{\widehat v_t}+\epsilon). $$
Products, squares, and divisions above are elementwise. AdamW decouples weight decay from the moment-adjusted step. For deep networks, important controls include learning-rate warmup/decay, gradient clipping, mixed-precision loss scaling where appropriate, and recording the effective batch size. None of these can substitute for a clean validation split.
4. C11 experiment: ridge regression by gradient descent
The following complete program trains $y\approx wx+b$ on a small dataset. It logs the mean squared error and incorporates ridge regularization. Bias is excluded from weight decay.
```c #include <stdio.h> #include <math.h>
int main(void) { const double x[] = {-2., -1., 0., 1., 2., 3.}; const double y[] = {-3., -1., 1., 3., 5., 7.}; const int n = (int)(sizeof(x) / sizeof(x[0])); const double eta = 0.04, lambda = 0.01; double w = 0., b = 0.; for (int t = 0; t < 6000; ++t) { double gw = 0., gb = 0.; for (int i = 0; i < n; ++i) { const double r = w*x[i] + b - y[i]; gw += 2.rx[i]/n; gb += 2.*r/n; } gw += 2.lambdaw; w -= eta*gw; b -= eta*gb; } double mse=0.; for (int i=0;i<n;++i) { double r=wx[i]+b-y[i];mse+=rr/n; } printf("w=%.6f b=%.6f mse=%.7f f(4)=%.6f\n",w,b,mse,4*w+b); return !(isfinite(mse) && mse<0.02); } ```
Compile: cc -std=c11 -O2 -Wall -Wextra ridge.c -lm -o ridge && ./ridge.
Because $\lambda>0$, the fitted slope is slightly below $2$; reducing $\lambda$ approaches the generating slope. This is a controlled toy problem and not a benchmark.
5. C++17 experiment: a differentiable two-layer XOR classifier
The XOR mapping $(0,0)\mapsto0$, $(1,0)\mapsto1$, $(0,1)\mapsto1$, $(1,1)\mapsto0$ is not linearly separable, so a single logistic neuron cannot represent it. A compact neural network with two nonlinear hidden units can. The code below performs explicit backpropagation over a 2–4–1 architecture, uses reproducible weights, and checks one parameter derivative numerically before training.
```cpp #include <array> #include <cmath> #include <iomanip> #include <iostream>
struct Point { std::array<double,2> x; double y; }; const std::array<Point,4> data{{{{0,0},0},{{0,1},1},{{1,0},1},{{1,1},0}}}; using Vec = std::array<double,17>; // W1(8), b1(4), W2(4), b2(1)
double sigmoid(double x) { if (x >= 0.) return 1. / (1. + std::exp(-x)); double e=std::exp(x); return e/(1.+e); } struct Result { double loss; Vec grad; }; Result evaluate(const Vec& p) { Result out{}; for (const auto& sample: data) { std::array<double,4> h{}; for (int j=0;j<4;++j) h[j]=std::tanh(p[2j]sample.x[0]+p[2j+1]sample.x[1]+p[8+j]); double logit=p[16]; for (int j=0;j<4;++j)logit+=p[12+j]*h[j]; double prob=sigmoid(logit); out.loss+=(std::max(0.,logit)-sample.y*logit+std::log1p(std::exp(-std::abs(logit))))/4.; double d=(prob-sample.y)/4.; for (int j=0;j<4;++j) { out.grad[12+j]+=d*h[j]; const double dh=dp[12+j](1.-h[j]*h[j]); out.grad[2j]+=dhsample.x[0]; out.grad[2j+1]+=dhsample.x[1]; out.grad[8+j]+=dh; } out.grad[16]+=d; } return out; } int main() { Vec p{{2.0,-1.5,-2.2,1.7,1.1,2.0,-1.4,-1.2, -0.5,-0.2,0.3,0.1,2.0,-1.9,1.3,0.8,0.0}}; auto g=evaluate(p); Vec plus=p,minus=p; const double eps=1e-5; plus[4]+=eps;minus[4]-=eps; double numeric=(evaluate(plus).loss-evaluate(minus).loss)/(2*eps); std::cout << "gradient discrepancy=" << std::abs(numeric-g.grad[4]) << '\n'; if (std::abs(numeric-g.grad[4])>1e-6)return 2; // Full-batch Adam with bias correction. Vec m{},v{};const double lr=0.045,beta1=0.9,beta2=0.999; for(int t=1;t<=6500;++t){ auto r=evaluate(p); for(int j=0;j<17;++j){ m[j]=beta1m[j]+(1-beta1)r.grad[j]; v[j]=beta2v[j]+(1-beta2)r.grad[j]*r.grad[j]; double mh=m[j]/(1-std::pow(beta1,t)); double vh=v[j]/(1-std::pow(beta2,t)); p[j]-=lr*mh/(std::sqrt(vh)+1e-8); } } std::cout << std::fixed << std::setprecision(5); for(const auto& ex:data){ double z=p[16]; for(int j=0;j<4;++j) z+=p[12+j]std::tanh(p[2j]ex.x[0]+p[2j+1]*ex.x[1]+p[8+j]); double prob=sigmoid(z); std::cout << ex.x[0] << "," << ex.x[1] << " -> " << prob << '\n'; if ((prob>0.5)!=bool(ex.y))return 3; } return 0; } ```
Compile: c++ -std=c++17 -O2 -Wall -Wextra xor.cpp -o xor && ./xor. The derivative check is a local implementation test. XOR is a tiny memorization task; success does not show out-of-distribution generalization.
6. Transformer pretraining: attention and objective functions
For a token representation matrix $X\in\mathbb R^{T\times d}$, the query, key and value projections are $Q=XW_Q$, $K=XW_K$, and $V=XW_V$. Scaled dot-product attention is
$$ \operatorname{Attention}(Q,K,V)=\operatorname{softmax}\!\left(\frac{QK^T}{\sqrt{d_k}}+M\right)V, $$
where a causal mask $M_{ij}=-\infty$ for $j>i$ prevents access to future tokens during autoregressive pretraining. In practice use stable softmax $s_i=\exp(a_i-m)/\sum_j\exp(a_j-m)$ with $m=\max_j a_j$. A basic attention matrix scales as $O(T^2)$ in sequence length; implementations such as memory-efficient tiled attention reduce the memory burden while preserving the mathematical operation up to numerical tolerances.
For next-token learning, let $x_1,\dots,x_T$ be token IDs. The causal language-model objective is
$$ L_{\mathrm{LM}}(\theta)=-\sum_{t=1}^{T-1}\log p_\theta(x_{t+1}\mid x_1,\ldots,x_t). $$
This is self-supervised because the text supplies its own labels. Predicting a token is not equivalent to verifying the truth of a claim. Overlap between training and evaluation documents can artificially inflate benchmark scores; deduplication and contamination audits matter.
6.1 Fine-tuning, low-rank adapters, and distillation
In supervised instruction tuning, one minimizes a task or response likelihood over curated prompt-response examples, sometimes masking prompt-token loss. In low-rank adaptation, a pretrained weight $W_0\in\mathbb R^{d_{out}\times d_{in}}$ is held fixed while
$$ W=W_0+\frac{\alpha}{r}BA,\qquad B\in\mathbb R^{d_{out}\times r},\quad A\in\mathbb R^{r\times d_{in}}. $$
The trainable adapter contains $r(d_{in}+d_{out})$ parameters rather than $d_{in}d_{out}$ for that layer. This can reduce optimizer-state requirements; it does not make the forward pass free. In response distillation a student may match a teacher distribution via $D_{KL}(p_T\Vert p_S)$, often combined with task loss; temperature scaling changes the smoothness and calibration of the soft targets.
6.2 Preferences, DPO, and objective alignment
Preference data contains pairs $(y^+,y^-)$ for a common prompt $x$. The common DPO formulation maximizes a logistic preference likelihood using a policy-versus-reference log-ratio:
$$ L_{\mathrm{DPO}}=-\mathbb E\left[\log\sigma\left(\beta\left(\log\frac{\pi_\theta(y^+|x)}{\pi_{ref}(y^+|x)}-\log\frac{\pi_\theta(y^-|x)}{\pi_{ref}(y^-|x)}\right)\right)\right]. $$
Preference learning optimizes the supplied preference signal; it is not a proof of safety or truthfulness. Effective evaluation includes adversarial tests, calibration where appropriate, and robustness across languages and domains.
7. Data governance, experiments, and reliable reporting
A defensible training pipeline documents data provenance, retention, deduplication, tokenization, train/validation/test boundaries and expected use rights. An RSS feed may provide title, date, link and short excerpt but not permission to republish a complete article or use its text for model training. A research archive and a training dataset therefore require different access rules.
When reporting a study, include: the baseline, random seeds, optimizer settings, number of updates, batch size, model size, device, wall-clock time, memory requirements, dataset licenses, and metrics with uncertainty estimates. If a confidence interval is claimed, report how it was computed (bootstrap, repeated runs or asymptotic approximation). For imbalanced classes, report precision, recall, confusion matrix and, when probabilities matter, calibration; accuracy alone is insufficient.
Use stable open-source implementations as baselines before asserting advantages. Ensure numerical comparisons use comparable hardware, preprocessing and stopping criteria. The two programs in this article are instructional implementations and do not benchmark production LLMs.
8. Recommended reproducible workflow
- Write the statistical hypothesis and define a measurable generalization metric before coding.
- Establish licensing and provenance for every data source; reject unauthorized training content.
- Split data without leakage and publish a baseline with an explicit random seed.
- Verify analytic gradients with finite differences on tiny deterministic problems.
- Scale batch sizes and model architecture only after validation curves and error categories justify them.
- Compare optimization choices, weight-decay behavior and stable-precision settings.
- Run repeated-seed and distribution-shift evaluations; publish failures, not only good examples.
- Export source-aware records in JSON/TXT/DOCX/PDF when legally appropriate; never confuse archival possession with training rights.
Conclusion
Contemporary AI training is not one secret algorithm. It combines statistical objectives, analytic differentiation, numerical optimization, representation learning and careful evaluation. The matrix derivations and executable C/C++ experiments show how these pieces fit together. Reliable scientific claims require reproducibility, proper data rights, quantitative baselines and an explicit distinction between an illustration and a validated result.
Selected primary references
- Rumelhart, Hinton and Williams, Learning representations by back-propagating errors (1986).
- Kingma and Ba, Adam: A Method for Stochastic Optimization (2014).
- Vaswani et al., Attention Is All You Need (2017).
- Loshchilov and Hutter, Decoupled Weight Decay Regularization (2017).
- Hu et al., LoRA: Low-Rank Adaptation of Large Language Models (2021).
- Rafailov et al., Direct Preference Optimization: Your Language Model is Secretly a Reward Model (2023).