Gradient Descent Convergence Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Strongly convex quadratic: f(x) = 0.5 * x^T A x - b^T x, with diagonal SPD A
5
// For diagonal A, L = max diag(A), mu = min diag(A). x* = A^{-1} b.
6
7
struct Quad {
8
    vector<double> diagA; // diagonal entries of A (positive)
9
    vector<double> b;     // vector b
10
11
    // Evaluate f(x)
12
    double f(const vector<double>& x) const {
13
        double v = 0.0;
14
        for (size_t i = 0; i < x.size(); ++i) {
15
            v += 0.5 * diagA[i] * x[i] * x[i] - b[i] * x[i];
16
        }
17
        return v;
18
    }
19
20
    // Gradient: A x - b (since A is diagonal)
21
    vector<double> grad(const vector<double>& x) const {
22
        vector<double> g(x.size());
23
        for (size_t i = 0; i < x.size(); ++i) g[i] = diagA[i] * x[i] - b[i];
24
        return g;
25
    }
26
27
    // Closed-form minimizer x* = A^{-1} b (element-wise since A is diagonal)
28
    vector<double> x_star() const {
29
        vector<double> xs(diagA.size());
30
        for (size_t i = 0; i < diagA.size(); ++i) xs[i] = b[i] / diagA[i];
31
        return xs;
32
    }
33
};
34
35
// Euclidean norm
36
double norm2(const vector<double>& x) {
37
    double s = 0.0;
38
    for (double v : x) s += v * v;
39
    return sqrt(s);
40
}
41
42
vector<double> sub(const vector<double>& a, const vector<double>& b) {
43
    vector<double> r(a.size());
44
    for (size_t i = 0; i < a.size(); ++i) r[i] = a[i] - b[i];
45
    return r;
46
}
47
48
int main() {
49
    ios::sync_with_stdio(false);
50
    cin.tie(nullptr);
51
52
    // Define a 5D problem with condition number kappa = L/mu = 10
53
    vector<double> diagA = {1.0, 2.0, 3.0, 5.0, 10.0}; // mu=1, L=10
54
    vector<double> b     = {1.0, -2.0, 0.5, 3.0, -1.0};
55
    Quad prob{diagA, b};
56
57
    const double L = *max_element(diagA.begin(), diagA.end());
58
    const double mu = *min_element(diagA.begin(), diagA.end());
59
    const double eta = 1.0 / L; // theoretical safe step size
60
61
    vector<double> x(5, 0.0);         // start at zero
62
    vector<double> xstar = prob.x_star();
63
    double dist0 = norm2(sub(x, xstar));
64
65
    cout << fixed << setprecision(6);
66
    cout << "# t\t||x_t - x*||\tBound (1 - mu/L)^t * ||x0 - x*||\t f(x_t) - f(x*)\n";
67
68
    // Precompute f(x*)
69
    double fxstar = prob.f(xstar);
70
71
    int T = 40;
72
    for (int t = 0; t <= T; ++t) {
73
        double bound = pow(1.0 - mu / L, t) * dist0; // linear convergence bound for eta=1/L
74
        double dist = norm2(sub(x, xstar));
75
        double gap = prob.f(x) - fxstar;
76
        cout << t << "\t" << dist << "\t" << bound << "\t" << gap << "\n";
77
78
        // Next iterate: x_{t+1} = x_t - eta * grad
79
        auto g = prob.grad(x);
80
        for (size_t i = 0; i < x.size(); ++i) x[i] -= eta * g[i];
81
    }
82
83
    return 0;
84
}
85

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Non-convex 1D function: f(x) = 0.25 * x^4 - 0.5 * x^2
5
// Gradient: f'(x) = x^3 - x; L-smooth on bounded sets; global non-convex.
6
7
int main() {
8
    ios::sync_with_stdio(false);
9
    cin.tie(nullptr);
10
11
    auto f = [](double x) { return 0.25 * x * x * x * x - 0.5 * x * x; };
12
    auto g = [](double x) { return x * x * x - x; };
13
14
    double x = 0.9;           // start near a non-convex region
15
    double eta = 0.2;         // modest step size to ensure descent
16
17
    cout << fixed << setprecision(6);
18
    cout << "# t\tx_t\tf(x_t)\t|grad|\tmin_{s<=t}|grad|^2 * t\n";
19
20
    double minGrad2 = numeric_limits<double>::infinity();
21
    int T = 200;
22
    for (int t = 1; t <= T; ++t) {
23
        double grad = g(x);
24
        minGrad2 = min(minGrad2, grad * grad);
25
        // Print a diagnostic capturing the O(1/t) stationarity trend
26
        cout << t << "\t" << x << "\t" << f(x) << "\t" << fabs(grad) << "\t" << (minGrad2 * t) << "\n";
27
        // Update
28
        x = x - eta * grad;
29
    }
30
31
    return 0;
32
}
33

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// 1D stochastic least squares: f(w) = (1/(2N)) sum_i (a_i*w - b_i)^2
5
// True minimizer: w* = (sum a_i b_i) / (sum a_i^2). We simulate SGD with eta_t = eta0 / sqrt(t).
6
7
int main() {
8
    ios::sync_with_stdio(false);
9
    cin.tie(nullptr);
10
11
    std::mt19937 rng(42);
12
    std::normal_distribution<double> noise(0.0, 0.5);
13
    std::uniform_real_distribution<double> a_dist(0.5, 2.0);
14
15
    const int N = 2000;
16
    const double w_true = 1.7;
17
18
    vector<double> a(N), b(N);
19
    for (int i = 0; i < N; ++i) {
20
        a[i] = a_dist(rng);
21
        b[i] = a[i] * w_true + noise(rng); // y = a w_true + noise
22
    }
23
24
    // Compute closed-form w* and f(w*) to quantify suboptimality
25
    double Sa2 = 0.0, Sab = 0.0;
26
    for (int i = 0; i < N; ++i) { Sa2 += a[i]*a[i]; Sab += a[i]*b[i]; }
27
    double w_star = Sab / Sa2;
28
    auto f = [&](double w){
29
        double s = 0.0;
30
        for (int i = 0; i < N; ++i) {
31
            double r = a[i]*w - b[i];
32
            s += 0.5 * r * r;
33
        }
34
        return s / N;
35
    };
36
    double f_star = f(w_star);
37
38
    // SGD with diminishing steps; sample one index per iteration
39
    std::uniform_int_distribution<int> uni(0, N-1);
40
41
    double w = 0.0;           // initialize
42
    double eta0 = 0.5;        // base step size (tuned)
43
44
    cout << fixed << setprecision(6);
45
    cout << "# t\tw_t\tf(w_t)-f*\t sqrt(t)*(f(w_t)-f*)\n";
46
47
    int T = 20000;
48
    for (int t = 1; t <= T; ++t) {
49
        int i = uni(rng);
50
        double grad = (a[i]*w - b[i]) * a[i]; // stochastic gradient from one sample
51
        double eta = eta0 / sqrt((double)t);
52
        w -= eta * grad;
53
54
        if (t % 1000 == 0) {
55
            double gap = max(0.0, f(w) - f_star);
56
            cout << t << "\t" << w << "\t" << gap << "\t" << (sqrt((double)t) * gap) << "\n";
57
        }
58
    }
59
60
    return 0;
61
}
62