Vectors & Vector Spaces

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <stdexcept>
5
#include <iomanip>
6
7
struct Vec {
8
    std::vector<double> a; // components
9
10
    explicit Vec(size_t n = 0, double val = 0.0) : a(n, val) {}
11
    explicit Vec(const std::vector<double>& v) : a(v) {}
12
13
    size_t size() const { return a.size(); }
14
15
    // Element access
16
    double& operator[](size_t i) { return a[i]; }
17
    const double& operator[](size_t i) const { return a[i]; }
18
19
    // Addition
20
    Vec operator+(const Vec& other) const {
21
        if (size() != other.size()) throw std::invalid_argument("dimension mismatch");
22
        Vec res(size());
23
        for (size_t i = 0; i < size(); ++i) res[i] = a[i] + other[i];
24
        return res;
25
    }
26
27
    // Subtraction
28
    Vec operator-(const Vec& other) const {
29
        if (size() != other.size()) throw std::invalid_argument("dimension mismatch");
30
        Vec res(size());
31
        for (size_t i = 0; i < size(); ++i) res[i] = a[i] - other[i];
32
        return res;
33
    }
34
35
    // Scalar multiplication
36
    Vec operator*(double s) const {
37
        Vec res(size());
38
        for (size_t i = 0; i < size(); ++i) res[i] = a[i] * s;
39
        return res;
40
    }
41
42
    // Dot product
43
    double dot(const Vec& other) const {
44
        if (size() != other.size()) throw std::invalid_argument("dimension mismatch");
45
        double sum = 0.0;
46
        for (size_t i = 0; i < size(); ++i) sum += a[i] * other[i];
47
        return sum;
48
    }
49
50
    // Euclidean norm
51
    double norm() const {
52
        return std::sqrt(this->dot(*this));
53
    }
54
55
    // Projection of this vector onto direction u: proj_u(this)
56
    Vec proj_onto(const Vec& u, double eps = 1e-12) const {
57
        double denom = u.dot(u);
58
        if (denom < eps) throw std::invalid_argument("cannot project onto near-zero vector");
59
        double coeff = this->dot(u) / denom;
60
        return u * coeff;
61
    }
62
};
63
64
std::ostream& operator<<(std::ostream& os, const Vec& v) {
65
    os << std::fixed << std::setprecision(4) << "[";
66
    for (size_t i = 0; i < v.size(); ++i) {
67
        os << v[i];
68
        if (i + 1 != v.size()) os << ", ";
69
    }
70
    os << "]";
71
    return os;
72
}
73
74
int main() {
75
    Vec v({1.0, 2.0, 3.0});
76
    Vec w({4.0, -1.0, 0.5});
77
78
    Vec sum = v + w;                 // addition
79
    Vec diff = v - w;                // subtraction
80
    Vec scaled = v * 3.0;            // scalar multiplication
81
    double dp = v.dot(w);            // dot product
82
    double nv = v.norm();            // norm
83
    Vec proj = v.proj_onto(w);       // projection of v onto w
84
85
    std::cout << "v        = " << v << "\n";
86
    std::cout << "w        = " << w << "\n";
87
    std::cout << "v + w    = " << sum << "\n";
88
    std::cout << "v - w    = " << diff << "\n";
89
    std::cout << "3 * v    = " << scaled << "\n";
90
    std::cout << "v·w      = " << dp << "\n";
91
    std::cout << "||v||    = " << nv << "\n";
92
    std::cout << "proj_w(v)= " << proj << "\n";
93
    return 0;
94
}
95

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Solve A x = b using Gaussian elimination with partial pivoting.
5
// A: d x k matrix (stored as vector<vector<double>> with d rows, k cols)
6
// b: size d
7
// Returns pair<has_solution, coefficients x of size k> (one solution if many)
8
pair<bool, vector<double>> solve_linear_system(vector<vector<double>> A, vector<double> b, double eps = 1e-10) {
9
    int d = (int)A.size();
10
    if (d == 0) return {true, {}};
11
    int k = (int)A[0].size();
12
    for (int i = 0; i < d; ++i) if ((int)A[i].size() != k) throw invalid_argument("inconsistent row sizes");
13
    if ((int)b.size() != d) throw invalid_argument("dimension mismatch");
14
15
    // Build augmented matrix [A | b]
16
    vector<vector<double>> M(d, vector<double>(k + 1));
17
    for (int i = 0; i < d; ++i) {
18
        for (int j = 0; j < k; ++j) M[i][j] = A[i][j];
19
        M[i][k] = b[i];
20
    }
21
22
    int row = 0, col = 0;
23
    vector<int> pivot_col; pivot_col.reserve(min(d, k));
24
25
    // Forward elimination to row-echelon form
26
    while (row < d && col < k) {
27
        // Choose pivot: row with largest |M[r][col]|
28
        int piv = row;
29
        for (int r = row + 1; r < d; ++r) if (fabs(M[r][col]) > fabs(M[piv][col])) piv = r;
30
        if (fabs(M[piv][col]) < eps) { // no pivot in this column
31
            ++col;
32
            continue;
33
        }
34
        swap(M[piv], M[row]);
35
        // Normalize pivot row
36
        double pivval = M[row][col];
37
        for (int j = col; j <= k; ++j) M[row][j] /= pivval;
38
        // Eliminate below
39
        for (int r = row + 1; r < d; ++r) {
40
            double factor = M[r][col];
41
            if (fabs(factor) < eps) continue;
42
            for (int j = col; j <= k; ++j) M[r][j] -= factor * M[row][j];
43
        }
44
        pivot_col.push_back(col);
45
        ++row; ++col;
46
    }
47
48
    // Check for inconsistency: a row of zeros in A but nonzero in b
49
    for (int r = row; r < d; ++r) {
50
        double sumA = 0.0;
51
        for (int j = 0; j < k; ++j) sumA += fabs(M[r][j]);
52
        if (sumA < eps && fabs(M[r][k]) > eps) return {false, {}}; // 0 = nonzero
53
    }
54
55
    // Back substitution to get one solution (set free vars = 0)
56
    vector<double> x(k, 0.0);
57
    for (int i = (int)pivot_col.size() - 1; i >= 0; --i) {
58
        int pc = pivot_col[i]; // column of pivot in row i
59
        double rhs = M[i][k];
60
        for (int j = pc + 1; j < k; ++j) rhs -= M[i][j] * x[j];
61
        x[pc] = rhs; // since pivot is 1 and others above are eliminated next
62
    }
63
64
    return {true, x};
65
}
66
67
// Convenience: check if v is in span of 'basis' columns, and retrieve coefficients
68
pair<bool, vector<double>> is_in_span(const vector<vector<double>>& basis, const vector<double>& v, double eps = 1e-10) {
69
    // basis is a list of vectors in R^d; form A with these as columns
70
    if (basis.empty()) return {all_of(v.begin(), v.end(), [](double x){ return fabs(x) < 1e-10; }), {}};
71
    int k = (int)basis.size();
72
    int d = (int)basis[0].size();
73
    for (const auto& b : basis) if ((int)b.size() != d) throw invalid_argument("basis vectors must have same dimension");
74
    if ((int)v.size() != d) throw invalid_argument("vector dimension mismatch");
75
76
    vector<vector<double>> A(d, vector<double>(k));
77
    for (int j = 0; j < k; ++j) for (int i = 0; i < d; ++i) A[i][j] = basis[j][i];
78
79
    return solve_linear_system(A, v, eps);
80
}
81
82
int main() {
83
    // Basis for the xy-plane in R^3
84
    vector<vector<double>> basis = { {1, 0, 0}, {0, 1, 0} };
85
    vector<double> v1 = {2, -3, 0};     // in span
86
    vector<double> v2 = {1, 1, 1};      // not in span
87
88
    auto r1 = is_in_span(basis, v1);
89
    auto r2 = is_in_span(basis, v2);
90
91
    cout << boolalpha;
92
    cout << "v1 in span? " << r1.first << "\n";
93
    if (r1.first) {
94
        cout << "coeffs for v1: [";
95
        for (size_t i = 0; i < r1.second.size(); ++i) cout << r1.second[i] << (i+1==r1.second.size()?"":", ");
96
        cout << "]\n";
97
    }
98
99
    cout << "v2 in span? " << r2.first << "\n";
100
    if (r2.first) {
101
        cout << "coeffs for v2: [";
102
        for (size_t i = 0; i < r2.second.size(); ++i) cout << r2.second[i] << (i+1==r2.second.size()?"":", ");
103
        cout << "]\n";
104
    } else {
105
        cout << "(no exact solution; z-component prevents membership)\n";
106
    }
107
    return 0;
108
}
109

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <algorithm>
5
#include <iomanip>
6
7
using namespace std;
8
9
static double dot(const vector<double>& a, const vector<double>& b) {
10
    if (a.size() != b.size()) throw invalid_argument("dimension mismatch");
11
    double s = 0.0; for (size_t i = 0; i < a.size(); ++i) s += a[i] * b[i];
12
    return s;
13
}
14
15
static double norm(const vector<double>& a) {
16
    return sqrt(dot(a, a));
17
}
18
19
static vector<double> add(const vector<double>& a, const vector<double>& b) {
20
    vector<double> r(a.size());
21
    for (size_t i = 0; i < a.size(); ++i) r[i] = a[i] + b[i];
22
    return r;
23
}
24
25
static vector<double> sub(const vector<double>& a, const vector<double>& b) {
26
    vector<double> r(a.size());
27
    for (size_t i = 0; i < a.size(); ++i) r[i] = a[i] - b[i];
28
    return r;
29
}
30
31
static vector<double> scal(double s, const vector<double>& a) {
32
    vector<double> r(a.size());
33
    for (size_t i = 0; i < a.size(); ++i) r[i] = s * a[i];
34
    return r;
35
}
36
37
// Modified Gram–Schmidt for numerical stability
38
vector<vector<double>> gram_schmidt_orthonormal(vector<vector<double>> vecs, double eps = 1e-12) {
39
    if (vecs.empty()) return {};
40
    size_t d = vecs[0].size();
41
    for (const auto& v : vecs) if (v.size() != d) throw invalid_argument("dimension mismatch");
42
43
    vector<vector<double>> U; U.reserve(vecs.size());
44
    for (auto v : vecs) {
45
        for (auto& u : U) {
46
            double c = dot(v, u);            // coefficient along u
47
            v = sub(v, scal(c, u));           // remove component along u
48
        }
49
        double nv = norm(v);
50
        if (nv > eps) {
51
            for (auto& x : v) x /= nv;       // normalize to unit length
52
            U.push_back(v);
53
        }
54
        // else: v was (nearly) dependent; skip it
55
    }
56
    return U; // orthonormal set spanning the same subspace
57
}
58
59
vector<double> project_onto_subspace(const vector<double>& v, const vector<vector<double>>& U) {
60
    if (U.empty()) return vector<double>(v.size(), 0.0);
61
    vector<double> p(v.size(), 0.0);
62
    for (const auto& u : U) {
63
        double c = dot(v, u);        // coefficient along orthonormal direction
64
        p = add(p, scal(c, u));
65
    }
66
    return p;
67
}
68
69
int main() {
70
    // Subspace spanned by two vectors in R^3
71
    vector<vector<double>> basis = { {1, 1, 0}, {1, 0, 1} };
72
    vector<double> v = {2, 1, 3};
73
74
    auto U = gram_schmidt_orthonormal(basis);
75
76
    cout << fixed << setprecision(6);
77
    cout << "Orthonormal basis U (rows):\n";
78
    for (auto& u : U) {
79
        cout << "[";
80
        for (size_t i = 0; i < u.size(); ++i) cout << u[i] << (i+1==u.size()?"":", ");
81
        cout << "]\n";
82
    }
83
84
    vector<double> p = project_onto_subspace(v, U);
85
    cout << "v = [" << v[0] << ", " << v[1] << ", " << v[2] << "]\n";
86
    cout << "Projection p = [" << p[0] << ", " << p[1] << ", " << p[2] << "]\n";
87
    cout << "Residual r = v - p = [" << v[0]-p[0] << ", " << v[1]-p[1] << ", " << v[2]-p[2] << "]\n";
88
    cout << "||r|| = " << norm(vector<double>{v[0]-p[0], v[1]-p[1], v[2]-p[2]}) << "\n";
89
90
    return 0;
91
}
92