Universal Approximation Theorem

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Simple 1D single-hidden-layer neural network with sigmoid activation.
5
// We train it by vanilla gradient descent to approximate f(x) = sin(2*pi*x) on [0,1].
6
7
struct Net1D {
8
    int H;                 // number of hidden units
9
    double c;              // output bias
10
    vector<double> w;      // input weights (size H)
11
    vector<double> b;      // hidden biases (size H)
12
    vector<double> a;      // output weights (size H)
13
14
    Net1D(int H_, unsigned seed=42): H(H_), c(0.0), w(H_), b(H_), a(H_) {
15
        mt19937 rng(seed);
16
        normal_distribution<double> nd(0.0, 0.5); // small random init
17
        for (int i = 0; i < H; ++i) {
18
            w[i] = nd(rng);
19
            b[i] = nd(rng);
20
            a[i] = nd(rng) * 0.1; // smaller output weights
21
        }
22
    }
23
24
    static inline double sigmoid(double z) {
25
        if (z >= 0) {
26
            double ez = exp(-z);
27
            return 1.0 / (1.0 + ez);
28
        } else {
29
            double ez = exp(z);
30
            return ez / (1.0 + ez);
31
        }
32
    }
33
34
    // Forward pass for a single scalar x
35
    double forward(double x, vector<double>* h_out=nullptr) const {
36
        vector<double> h(H);
37
        for (int i = 0; i < H; ++i) {
38
            h[i] = sigmoid(w[i]*x + b[i]);
39
        }
40
        double y = c;
41
        for (int i = 0; i < H; ++i) y += a[i]*h[i];
42
        if (h_out) *h_out = move(h);
43
        return y;
44
    }
45
};
46
47
static inline double target(double x) {
48
    const double PI = acos(-1.0);
49
    return sin(2.0 * PI * x);
50
}
51
52
int main() {
53
    ios::sync_with_stdio(false);
54
    cin.tie(nullptr);
55
56
    int H = 50;                 // hidden units
57
    int N = 256;                // training samples
58
    int epochs = 2000;          // training epochs
59
    double lr = 0.05;           // learning rate
60
61
    // Generate training data uniformly on [0,1]
62
    vector<double> xs(N), ys(N);
63
    for (int i = 0; i < N; ++i) {
64
        xs[i] = (i + 0.5) / N;
65
        ys[i] = target(xs[i]);
66
    }
67
68
    Net1D net(H, 123);
69
70
    // Training loop: mean squared error with full-batch gradient descent
71
    for (int epoch = 0; epoch < epochs; ++epoch) {
72
        // Accumulate gradients
73
        double gc = 0.0;              // grad wrt c
74
        vector<double> ga(H, 0.0);    // grad wrt a
75
        vector<double> gw(H, 0.0);    // grad wrt w
76
        vector<double> gb(H, 0.0);    // grad wrt b
77
78
        double mse = 0.0;
79
        for (int i = 0; i < N; ++i) {
80
            vector<double> h;
81
            double yhat = net.forward(xs[i], &h);
82
            double err = yhat - ys[i];
83
            mse += err * err;
84
85
            // dL/dc = 2/N * err
86
            gc += 2.0 * err / N;
87
            // dL/da_j = 2/N * err * h_j
88
            for (int j = 0; j < net.H; ++j) ga[j] += 2.0 * err * h[j] / N;
89
            // For sigmoid, dh_j/dz = h_j*(1-h_j); z = w_j*x + b_j
90
            for (int j = 0; j < net.H; ++j) {
91
                double dh_dz = h[j] * (1.0 - h[j]);
92
                double dz_dw = xs[i];
93
                double dz_db = 1.0;
94
                double dL_dz = 2.0 * err * net.a[j] * dh_dz / N;
95
                gw[j] += dL_dz * dz_dw;
96
                gb[j] += dL_dz * dz_db;
97
            }
98
        }
99
        mse /= N;
100
101
        // Gradient step
102
        net.c -= lr * gc;
103
        for (int j = 0; j < net.H; ++j) {
104
            net.a[j] -= lr * ga[j];
105
            net.w[j] -= lr * gw[j];
106
            net.b[j] -= lr * gb[j];
107
        }
108
109
        if (epoch % 200 == 0) {
110
            cout << "Epoch " << epoch << ", MSE = " << mse << "\n";
111
        }
112
    }
113
114
    // Report a few predictions
115
    cout << fixed << setprecision(6);
116
    for (double x : {0.0, 0.125, 0.25, 0.5, 0.75, 1.0}) {
117
        double yhat = net.forward(x);
118
        cout << "x=" << setw(6) << x << "  f(x)=" << setw(9) << target(x)
119
             << "  yhat=" << setw(10) << yhat << "\n";
120
    }
121
122
    return 0;
123
}
124

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Approximate a 2D function using a single hidden layer with fixed random ReLU features
5
// and solve the output weights via ridge regression: (Z^T Z + λI)a = Z^T y.
6
7
struct RandomFeatures2D {
8
    int H; // number of hidden units
9
    vector<array<double,2>> W; // hidden weights w_i in R^2
10
    vector<double> b;          // hidden biases
11
    RandomFeatures2D(int H_, unsigned seed=42): H(H_), W(H_), b(H_) {
12
        mt19937 rng(seed);
13
        normal_distribution<double> nd(0.0, 2.0); // somewhat wide to cover [0,1]^2
14
        uniform_real_distribution<double> ub(-1.0, 1.0);
15
        for (int i = 0; i < H; ++i) {
16
            W[i] = {nd(rng), nd(rng)};
17
            b[i] = ub(rng);
18
        }
19
    }
20
    static inline double relu(double z) { return z > 0 ? z : 0; }
21
    // Compute hidden features for a single (x,y)
22
    vector<double> phi(double x, double y) const {
23
        vector<double> h(H);
24
        for (int i = 0; i < H; ++i) {
25
            double z = W[i][0]*x + W[i][1]*y + b[i];
26
            h[i] = relu(z);
27
        }
28
        return h;
29
    }
30
};
31
32
// Solve (A)x = b by Gaussian elimination with partial pivoting (small H only)
33
bool gaussian_solve(vector<vector<double>>& A, vector<double>& b, vector<double>& x) {
34
    int n = (int)A.size();
35
    x.assign(n, 0.0);
36
    // Augment A with b
37
    for (int i = 0; i < n; ++i) A[i].push_back(b[i]);
38
    // Elimination
39
    for (int col = 0; col < n; ++col) {
40
        // Pivot
41
        int piv = col;
42
        for (int r = col+1; r < n; ++r)
43
            if (fabs(A[r][col]) > fabs(A[piv][col])) piv = r;
44
        if (fabs(A[piv][col]) < 1e-12) return false; // singular
45
        swap(A[piv], A[col]);
46
        // Normalize row
47
        double div = A[col][col];
48
        for (int j = col; j <= n; ++j) A[col][j] /= div;
49
        // Eliminate others
50
        for (int r = 0; r < n; ++r) if (r != col) {
51
            double factor = A[r][col];
52
            for (int j = col; j <= n; ++j) A[r][j] -= factor * A[col][j];
53
        }
54
    }
55
    for (int i = 0; i < n; ++i) x[i] = A[i][n];
56
    return true;
57
}
58
59
static inline double ftarget(double x, double y) {
60
    // A smooth 2D bump centered at (0.5, 0.5)
61
    double dx = x - 0.5, dy = y - 0.5;
62
    return exp(-(dx*dx + dy*dy) / 0.02);
63
}
64
65
int main() {
66
    ios::sync_with_stdio(false);
67
    cin.tie(nullptr);
68
69
    int H = 200;          // random ReLU features
70
    double lambda = 1e-3; // ridge regularization
71
72
    // Build training grid on [0,1]^2
73
    int G = 20; // 20x20 grid => N=400
74
    int N = G * G;
75
    vector<array<double,2>> X(N);
76
    vector<double> y(N);
77
    int idx = 0;
78
    for (int i = 0; i < G; ++i) {
79
        for (int j = 0; j < G; ++j) {
80
            double x = (i + 0.5) / G;
81
            double yy = (j + 0.5) / G;
82
            X[idx] = {x, yy};
83
            y[idx] = ftarget(x, yy);
84
            ++idx;
85
        }
86
    }
87
88
    RandomFeatures2D rf(H, 7);
89
90
    // Compute feature matrix Z (N x H)
91
    vector<vector<double>> Z(N, vector<double>(H));
92
    for (int i = 0; i < N; ++i) {
93
        auto h = rf.phi(X[i][0], X[i][1]);
94
        for (int j = 0; j < H; ++j) Z[i][j] = h[j];
95
    }
96
97
    // Form normal equations: (Z^T Z + λI)a = Z^T y
98
    vector<vector<double>> A(H, vector<double>(H, 0.0));
99
    vector<double> By(H, 0.0);
100
    for (int i = 0; i < H; ++i) {
101
        for (int j = 0; j < H; ++j) {
102
            double s = 0.0;
103
            for (int k = 0; k < N; ++k) s += Z[k][i] * Z[k][j];
104
            A[i][j] = s + (i == j ? lambda : 0.0);
105
        }
106
        double t = 0.0;
107
        for (int k = 0; k < N; ++k) t += Z[k][i] * y[k];
108
        By[i] = t;
109
    }
110
111
    // Solve for output weights a
112
    vector<double> a;
113
    bool ok = gaussian_solve(A, By, a);
114
    if (!ok) {
115
        cerr << "Failed to solve normal equations (possibly ill-conditioned).\n";
116
        return 0;
117
    }
118
119
    // Evaluate training MSE
120
    double mse = 0.0;
121
    for (int i = 0; i < N; ++i) {
122
        double pred = 0.0;
123
        for (int j = 0; j < H; ++j) pred += Z[i][j] * a[j];
124
        double err = pred - y[i];
125
        mse += err * err;
126
    }
127
    mse /= N;
128
    cout << fixed << setprecision(6);
129
    cout << "Training MSE: " << mse << "\n";
130
131
    // Test a few points
132
    for (auto P : vector<array<double,2>>{{0.1,0.1},{0.5,0.5},{0.8,0.2}}) {
133
        auto h = rf.phi(P[0], P[1]);
134
        double pred = 0.0;
135
        for (int j = 0; j < H; ++j) pred += h[j] * a[j];
136
        cout << "(x,y)=(" << P[0] << "," << P[1] << ")  f=" << ftarget(P[0],P[1])
137
             << "  yhat=" << pred << "\n";
138
    }
139
140
    return 0;
141
}
142

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Approximate a 1D continuous function on [0,1] by linear splines and
5
// represent the spline exactly as a ReLU network: g(x) = c0 + c1*x + sum a_j * ReLU(x - t_j).
6
// This demonstrates a constructive instance of the UAT in 1D.
7
8
static inline double relu(double z) { return z > 0 ? z : 0; }
9
10
// Build spline coefficients from sampled points (t_0=0 < t_1 < ... < t_K=1)
11
struct ReLUSpline {
12
    vector<double> t;   // breakpoints t_j, size K+1
13
    vector<double> a;   // hinge weights a_j, size K-1 (for j=1..K-1)
14
    double c0;          // intercept
15
    double c1;          // slope of first segment
16
17
    // Given samples (t_j, f_j), produce exact PWL interpolant and convert to ReLU form
18
    ReLUSpline(const vector<double>& T, const vector<double>& F) {
19
        int K = (int)T.size() - 1; // number of segments
20
        t = T;
21
        // Compute segment slopes s_j on [t_j, t_{j+1}]
22
        vector<double> s(K);
23
        for (int j = 0; j < K; ++j) {
24
            s[j] = (F[j+1] - F[j]) / (T[j+1] - T[j]);
25
        }
26
        c0 = F[0];
27
        c1 = s[0];
28
        // a_j are slope jumps: a_j = s_j - s_{j-1} for j=1..K-1
29
        a.assign(K-1, 0.0);
30
        for (int j = 1; j < K; ++j) a[j-1] = s[j] - s[j-1];
31
    }
32
33
    double eval(double x) const {
34
        // g(x) = c0 + c1*x + sum_{j=1}^{K-1} a_j * ReLU(x - t_j)
35
        double y = c0 + c1 * x;
36
        for (int j = 0; j < (int)a.size(); ++j) y += a[j] * relu(x - t[j+1]);
37
        return y;
38
    }
39
};
40
41
static inline double ftrue(double x) {
42
    // target function to approximate
43
    const double PI = acos(-1.0);
44
    return 0.5 * sin(2.0 * PI * x) + 0.3 * x * (1.0 - x);
45
}
46
47
int main() {
48
    ios::sync_with_stdio(false);
49
    cin.tie(nullptr);
50
51
    // Choose K segments for the spline
52
    int K = 16; // more segments => better uniform approximation
53
    vector<double> T(K+1), F(K+1);
54
    for (int j = 0; j <= K; ++j) {
55
        T[j] = (double)j / K;
56
        F[j] = ftrue(T[j]);
57
    }
58
59
    ReLUSpline spline(T, F);
60
61
    // Evaluate error on a fine grid
62
    int M = 1000;
63
    double maxerr = 0.0;
64
    for (int i = 0; i <= M; ++i) {
65
        double x = (double)i / M;
66
        double err = fabs(spline.eval(x) - ftrue(x));
67
        if (err > maxerr) maxerr = err;
68
    }
69
70
    cout << fixed << setprecision(6);
71
    cout << "Uniform error with K=" << K << " segments: " << maxerr << "\n";
72
73
    // Show a few sample evaluations
74
    for (double x : {0.0, 0.1, 0.25, 0.5, 0.75, 1.0}) {
75
        cout << "x=" << x << "  f(x)=" << ftrue(x) << "  spline(x)=" << spline.eval(x) << "\n";
76
    }
77
78
    return 0;
79
}
80