Hilbert Spaces

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Compute inner product (dot product) in R^n
5
double inner_product_vec(const vector<double>& x, const vector<double>& y) {
6
    assert(x.size() == y.size());
7
    // Kahan compensated summation to reduce floating-point error
8
    double sum = 0.0, c = 0.0;
9
    for (size_t i = 0; i < x.size(); ++i) {
10
        double prod = x[i] * y[i] - c;
11
        double t = sum + prod;
12
        c = (t - sum) - prod;
13
        sum = t;
14
    }
15
    return sum;
16
}
17
18
// Euclidean norm induced by the inner product
19
double norm_vec(const vector<double>& x) {
20
    return sqrt(max(0.0, inner_product_vec(x, x)));
21
}
22
23
// Project y onto the (closed) subspace span{x}. If x is the zero vector, returns zero vector.
24
vector<double> project_onto_vector(const vector<double>& y, const vector<double>& x) {
25
    vector<double> proj(y.size(), 0.0);
26
    double denom = inner_product_vec(x, x);
27
    if (denom == 0.0) return proj; // projection onto {0}
28
    double alpha = inner_product_vec(y, x) / denom; // coefficient ⟨y,x⟩ / ⟨x,x⟩
29
    for (size_t i = 0; i < y.size(); ++i) proj[i] = alpha * x[i];
30
    return proj;
31
}
32
33
int main() {
34
    // Example vectors in R^5
35
    vector<double> x = {1, 2, 0, -1, 3};
36
    vector<double> y = {4, 0, 1, 2, -2};
37
38
    cout << fixed << setprecision(6);
39
    cout << "<x,y> = " << inner_product_vec(x, y) << "\n";
40
    cout << "||x|| = " << norm_vec(x) << ", ||y|| = " << norm_vec(y) << "\n";
41
42
    vector<double> p = project_onto_vector(y, x);
43
    // Residual r = y - p is orthogonal to x in exact arithmetic
44
    vector<double> r(y.size());
45
    for (size_t i = 0; i < y.size(); ++i) r[i] = y[i] - p[i];
46
47
    cout << "Projection of y onto span{x}: [";
48
    for (size_t i = 0; i < p.size(); ++i) cout << p[i] << (i+1==p.size()? "]\n" : ", ");
49
50
    cout << "Orthogonality check <x, y - P x> = " << inner_product_vec(x, r) << " (should be ~0)\n";
51
    return 0;
52
}
53

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Inner product and norm as before
5
double inner_product_vec(const vector<double>& x, const vector<double>& y) {
6
    assert(x.size() == y.size());
7
    double sum = 0.0, c = 0.0;
8
    for (size_t i = 0; i < x.size(); ++i) {
9
        double prod = x[i] * y[i] - c;
10
        double t = sum + prod;
11
        c = (t - sum) - prod;
12
        sum = t;
13
    }
14
    return sum;
15
}
16
17
double norm_vec(const vector<double>& x) { return sqrt(max(0.0, inner_product_vec(x, x))); }
18
19
// Modified Gram-Schmidt: orthonormalize columns in V (vectors of dimension n)
20
vector<vector<double>> modified_gram_schmidt(const vector<vector<double>>& V) {
21
    size_t m = V.size();
22
    if (m == 0) return {};
23
    size_t n = V[0].size();
24
    for (size_t i = 1; i < m; ++i) assert(V[i].size() == n);
25
26
    vector<vector<double>> Q; // orthonormal vectors
27
    Q.reserve(m);
28
29
    for (size_t j = 0; j < m; ++j) {
30
        vector<double> v = V[j]; // copy to modify
31
        // subtract projections onto previous q_i
32
        for (size_t i = 0; i < Q.size(); ++i) {
33
            double rij = inner_product_vec(Q[i], v);
34
            for (size_t k = 0; k < n; ++k) v[k] -= rij * Q[i][k];
35
        }
36
        double normv = norm_vec(v);
37
        if (normv > 0) {
38
            for (size_t k = 0; k < n; ++k) v[k] /= normv; // normalize
39
            Q.push_back(move(v));
40
        } else {
41
            // dependent vector; skip
42
        }
43
    }
44
    return Q;
45
}
46
47
// Project vector b onto span(Q) where Q is orthonormal
48
vector<double> project_onto_subspace(const vector<double>& b, const vector<vector<double>>& Q) {
49
    size_t n = b.size();
50
    vector<double> p(n, 0.0);
51
    for (const auto& q : Q) {
52
        double coeff = inner_product_vec(b, q); // ⟨b, q⟩
53
        for (size_t i = 0; i < n; ++i) p[i] += coeff * q[i];
54
    }
55
    return p;
56
}
57
58
int main() {
59
    // Columns (as vectors) of a design matrix in R^4
60
    vector<vector<double>> V = {
61
        {1, 1, 1, 1},     // constant
62
        {-1, 0, 1, 2},    // linear-ish
63
        {1, 0, 1, 4}      // quadratic-ish
64
    };
65
    // Orthonormalize columns
66
    auto Q = modified_gram_schmidt(V);
67
68
    // Vector to approximate (e.g., measurements)
69
    vector<double> b = {0.9, 1.0, 2.2, 2.8};
70
71
    auto p = project_onto_subspace(b, Q); // least-squares approximation
72
73
    cout << fixed << setprecision(6);
74
    cout << "Orthonormal basis size: " << Q.size() << "\n";
75
    // Check orthonormality
76
    for (size_t i = 0; i < Q.size(); ++i) {
77
        for (size_t j = i; j < Q.size(); ++j) {
78
            cout << "<q" << i << ", q" << j << "> = " << inner_product_vec(Q[i], Q[j]) << "\n";
79
        }
80
    }
81
82
    // Report projection and residual norm
83
    double res_norm = 0.0;
84
    cout << "Projection p = [";
85
    for (size_t i = 0; i < p.size(); ++i) {
86
        cout << p[i] << (i + 1 == p.size() ? "]\n" : ", ");
87
    }
88
    vector<double> r(b.size());
89
    for (size_t i = 0; i < b.size(); ++i) r[i] = b[i] - p[i];
90
    cout << "Residual norm ||b - p|| = " << norm_vec(r) << "\n";
91
    return 0;
92
}
93

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Approximate the inner product <f,g> = (1/pi) * ∫_0^{2π} f(x) g(x) dx using the trapezoidal rule.
5
// For periodic smooth functions, this rule is highly accurate on uniform grids.
6
const double PI = acos(-1.0);
7
8
double inner_product_periodic(function<double(double)> f,
9
                              function<double(double)> g,
10
                              int N) {
11
    // N uniform samples on [0, 2π]
12
    double h = 2.0 * PI / N;
13
    auto w = [&](int i){ return (i==0 || i==N) ? 0.5 : 1.0; }; // trapezoid weights at endpoints
14
15
    double sum = 0.0;
16
    for (int i = 0; i <= N; ++i) {
17
        double x = i * h;
18
        double fx = f(x);
19
        double gx = g(x);
20
        sum += w(i) * fx * gx;
21
    }
22
    return (1.0 / PI) * h * sum; // (1/π) * integral approximation
23
}
24
25
// Orthonormal trigonometric system on [0, 2π] under <f,g> = (1/π)∫ f g : {1/√2, cos(kx), sin(kx)}
26
struct FourierCoeffs {
27
    double a0_over_sqrt2; // coefficient for 1/√2
28
    vector<double> a;     // cos coefficients
29
    vector<double> b;     // sin coefficients
30
};
31
32
FourierCoeffs project_fourier(function<double(double)> f, int K, int Nsamples) {
33
    FourierCoeffs c; c.a.resize(K+1); c.b.resize(K+1);
34
    // constant basis e0 = 1/√2
35
    c.a0_over_sqrt2 = inner_product_periodic(f, [](double){ return 1.0 / sqrt(2.0); }, Nsamples);
36
    for (int k = 1; k <= K; ++k) {
37
        auto cosk = [k](double x){ return cos(k * x); };
38
        auto sink = [k](double x){ return sin(k * x); };
39
        c.a[k] = inner_product_periodic(f, cosk, Nsamples);
40
        c.b[k] = inner_product_periodic(f, sink, Nsamples);
41
    }
42
    return c;
43
}
44
45
// Evaluate truncated Fourier series S_K at x
46
double fourier_partial_sum(const FourierCoeffs& c, double x) {
47
    int K = (int)c.a.size() - 1;
48
    double s = c.a0_over_sqrt2 / sqrt(2.0); // a0 * (1/√2)
49
    for (int k = 1; k <= K; ++k) s += c.a[k] * cos(k*x) + c.b[k] * sin(k*x);
50
    return s;
51
}
52
53
int main() {
54
    cout << fixed << setprecision(6);
55
56
    // Target function in L^2(0,2π): f(x) = x (not periodic in value, but square-integrable; expect Gibbs near endpoints)
57
    auto f = [](double x){ return x; };
58
59
    int K = 10;        // number of Fourier modes
60
    int N = 4096;      // quadrature samples (even, large improves accuracy)
61
62
    auto coeffs = project_fourier(f, K, N);
63
64
    // Report a few coefficients
65
    cout << "a0/√2 = " << coeffs.a0_over_sqrt2 << "\n";
66
    for (int k = 1; k <= min(K, 5); ++k) {
67
        cout << "a[" << k << "] = " << coeffs.a[k] << ", b[" << k << "] = " << coeffs.b[k] << "\n";
68
    }
69
70
    // Compute L^2-like error on grid: sqrt((1/π)∫ |f - S_K|^2)
71
    double h = 2.0 * PI / N;
72
    double sum = 0.0;
73
    for (int i = 0; i <= N; ++i) {
74
        double x = i * h;
75
        double fx = f(x);
76
        double sk = fourier_partial_sum(coeffs, x);
77
        double w = (i==0 || i==N) ? 0.5 : 1.0; // trapezoid weight
78
        sum += w * (fx - sk) * (fx - sk);
79
    }
80
    double l2err = sqrt((1.0 / PI) * h * sum);
81
    cout << "Approximate L2 error with K=" << K << " is " << l2err << "\n";
82
83
    // Show reconstruction at a few points
84
    for (double x : {0.0, PI/2, PI, 3*PI/2, 2*PI}) {
85
        cout << "x=" << x << ": f(x)=" << f(x) << ", S_K(x)=" << fourier_partial_sum(coeffs, x) << "\n";
86
    }
87
88
    return 0;
89
}
90