Universal Approximation Theorems

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// A simple 1D -> 1D MLP with one hidden layer and sigmoid activation
5
struct MLP1D {
6
    int H;                 // number of hidden units
7
    vector<double> W;      // hidden weights (size H), since input is scalar
8
    vector<double> b;      // hidden biases (size H)
9
    vector<double> a;      // output weights (size H)
10
    double beta;           // output bias
11
    mt19937 rng;
12
13
    MLP1D(int H_, uint32_t seed = 42) : H(H_), W(H_), b(H_), a(H_), beta(0.0), rng(seed) {
14
        normal_distribution<double> nd(0.0, 0.5); // small random init
15
        for (int i = 0; i < H; ++i) {
16
            W[i] = nd(rng);
17
            b[i] = nd(rng);
18
            a[i] = nd(rng);
19
        }
20
        beta = 0.0;
21
    }
22
23
    static inline double sigmoid(double x) {
24
        // Numerically stable sigmoid
25
        if (x >= 0) {
26
            double z = exp(-x);
27
            return 1.0 / (1.0 + z);
28
        } else {
29
            double z = exp(x);
30
            return z / (1.0 + z);
31
        }
32
    }
33
34
    // Forward pass for a single scalar x; returns y_hat and stores hidden activations if provided
35
    double forward(double x, vector<double>* h_ptr = nullptr) const {
36
        vector<double> hlocal;
37
        vector<double>& h = h_ptr ? *h_ptr : hlocal;
38
        if (!h_ptr) h.assign(H, 0.0);
39
        double y = beta;
40
        for (int i = 0; i < H; ++i) {
41
            double z = W[i] * x + b[i];
42
            double s = sigmoid(z);
43
            if (!h_ptr) h[i] = s; else h[i] = s;
44
            y += a[i] * s;
45
        }
46
        return y;
47
    }
48
49
    // One SGD step on a single (x, y) with learning rate eta
50
    void sgd_step(double x, double y, double eta) {
51
        vector<double> h(H);
52
        double y_hat = forward(x, &h);
53
        double diff = y_hat - y; // dL/dy_hat for MSE with factor 1
54
55
        // Gradients
56
        double dbeta = diff; // dL/dbeta = diff
57
        vector<double> da(H), dW(H), dbv(H);
58
        for (int i = 0; i < H; ++i) {
59
            da[i] = diff * h[i];
60
            // backprop through sigmoid: h = sigmoid(z), dh/dz = h*(1-h)
61
            double dz = diff * a[i] * h[i] * (1.0 - h[i]);
62
            dW[i] = dz * x; // since z = W[i]*x + b[i]
63
            dbv[i] = dz;
64
        }
65
66
        // Parameter update
67
        beta -= eta * dbeta;
68
        for (int i = 0; i < H; ++i) {
69
            a[i]   -= eta * da[i];
70
            W[i]   -= eta * dW[i];
71
            b[i]   -= eta * dbv[i];
72
        }
73
    }
74
};
75
76
int main() {
77
    ios::sync_with_stdio(false);
78
    cin.tie(nullptr);
79
80
    // Target function f(x) = sin(2*pi*x)
81
    auto f = [](double x){ return sin(2.0 * M_PI * x); };
82
83
    // Generate training data on [0,1]
84
    const int N = 200;
85
    vector<double> X(N), Y(N);
86
    std::mt19937 rng(123);
87
    std::uniform_real_distribution<double> unif(0.0, 1.0);
88
    for (int i = 0; i < N; ++i) {
89
        X[i] = unif(rng);
90
        Y[i] = f(X[i]);
91
    }
92
93
    // Initialize model
94
    int H = 32;           // hidden width; increase for tighter fits
95
    double eta = 0.1;     // learning rate
96
    int epochs = 2000;    // number of passes over data
97
98
    MLP1D net(H, 999);
99
100
    // Training loop: shuffle each epoch for SGD
101
    vector<int> idx(N); iota(idx.begin(), idx.end(), 0);
102
    for (int e = 0; e < epochs; ++e) {
103
        shuffle(idx.begin(), idx.end(), rng);
104
        for (int id : idx) {
105
            net.sgd_step(X[id], Y[id], eta);
106
        }
107
        // Monitor every 200 epochs
108
        if ((e+1) % 200 == 0) {
109
            // Compute MSE on a dense grid
110
            double mse = 0.0; int M = 1000;
111
            for (int i = 0; i < M; ++i) {
112
                double x = (i + 0.5) / M;
113
                double y_hat = net.forward(x);
114
                double d = y_hat - f(x);
115
                mse += d * d;
116
            }
117
            mse /= M;
118
            cerr << "Epoch " << (e+1) << ": MSE(grid) = " << mse << "\n";
119
        }
120
    }
121
122
    // Print a few sample approximations
123
    cout << fixed << setprecision(6);
124
    cout << "x\tf(x)\ty_hat(x)\n";
125
    for (int i = 0; i <= 10; ++i) {
126
        double x = i / 10.0;
127
        cout << x << "\t" << f(x) << "\t" << net.forward(x) << "\n";
128
    }
129
130
    return 0;
131
}
132

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Given knots (x0 < x1 < ... < xK) and values y_i, build
5
// g(x) = y0 + m0*x + sum_{i=1}^{K-1} (m_i - m_{i-1}) * ReLU(x - x_i),
6
// where m_i = (y_{i+1} - y_i) / (x_{i+1} - x_i). This equals the linear interpolant.
7
8
struct ReLUInterp {
9
    vector<double> t;   // thresholds x_i, i = 1..K-1
10
    vector<double> c;   // coefficients (m_i - m_{i-1}), size K-1
11
    double y0;          // base value
12
    double m0;          // initial slope
13
14
    static inline double relu(double x) { return x > 0 ? x : 0; }
15
16
    // Build from knots x[0..K], y[0..K], with strictly increasing x
17
    ReLUInterp(const vector<double>& x, const vector<double>& y) {
18
        int K = (int)x.size() - 1;
19
        if ((int)y.size() != K + 1 || K < 1) {
20
            throw runtime_error("Invalid input sizes");
21
        }
22
        for (int i = 0; i < K; ++i) {
23
            if (!(x[i+1] > x[i])) throw runtime_error("x must be strictly increasing");
24
        }
25
        vector<double> m(K); // slopes on each interval [x_i, x_{i+1}]
26
        for (int i = 0; i < K; ++i) {
27
            m[i] = (y[i+1] - y[i]) / (x[i+1] - x[i]);
28
        }
29
        y0 = y[0];
30
        m0 = m[0];
31
        t.clear(); c.clear();
32
        for (int i = 1; i < K; ++i) {
33
            t.push_back(x[i]);
34
            c.push_back(m[i] - m[i-1]);
35
        }
36
    }
37
38
    // Evaluate the network at x
39
    double operator()(double xval) const {
40
        double g = y0 + m0 * xval;
41
        for (size_t i = 0; i < t.size(); ++i) {
42
            g += c[i] * relu(xval - t[i]);
43
        }
44
        return g;
45
    }
46
};
47
48
int main() {
49
    ios::sync_with_stdio(false);
50
    cin.tie(nullptr);
51
52
    // Define knots and values of a smooth function, e.g., f(x) = exp(-x) * cos(3x)
53
    auto f = [](double x){ return exp(-x) * cos(3.0 * x); };
54
55
    int K = 8; // number of segments
56
    vector<double> x(K+1), y(K+1);
57
    double xmin = 0.0, xmax = 2.0;
58
    for (int i = 0; i <= K; ++i) {
59
        x[i] = xmin + (xmax - xmin) * i / K;
60
        y[i] = f(x[i]);
61
    }
62
63
    // Build the ReLU network representing the linear interpolant
64
    ReLUInterp net(x, y);
65
66
    // Verify exactness at knots and print max error on a dense grid
67
    cout << fixed << setprecision(6);
68
    cout << "Knot checks (x, f(x), g(x))\n";
69
    for (int i = 0; i <= K; ++i) {
70
        double gx = net(x[i]);
71
        cout << x[i] << "\t" << y[i] << "\t" << gx << "\n";
72
    }
73
74
    // Evaluate error on a dense grid
75
    double max_err = 0.0; int M = 2000;
76
    for (int i = 0; i <= M; ++i) {
77
        double xv = xmin + (xmax - xmin) * i / M;
78
        double err = fabs(net(xv) - f(xv));
79
        if (err > max_err) max_err = err;
80
    }
81
    cerr << "Max error vs. true f on grid (should decrease as K increases): " << max_err << "\n";
82
83
    return 0;
84
}
85