Partial Derivatives

1
#include <cmath>
2
#include <functional>
3
#include <iostream>
4
#include <limits>
5
#include <vector>
6
7
// Compute the i-th partial derivative of f at point x using central differences.
8
double partial_derivative(const std::function<double(const std::vector<double>&)>& f,
9
                          std::vector<double> x,
10
                          std::size_t i) {
11
    const double xi = x[i];
12
    const double eps = std::numeric_limits<double>::epsilon();
13
    // Heuristic step size: h ~ sqrt(eps)*max(1, |xi|)
14
    const double h = std::sqrt(eps) * std::max(1.0, std::abs(xi));
15
16
    x[i] = xi + h;
17
    const double fp = f(x);
18
    x[i] = xi - h;
19
    const double fm = f(x);
20
21
    // Restore (optional)
22
    x[i] = xi;
23
24
    return (fp - fm) / (2.0 * h);
25
}
26
27
int main() {
28
    // Example scalar function: f(x, y, z) = x^2*y + sin(y) + x*exp(z)
29
    auto f = [](const std::vector<double>& v) -> double {
30
        double x = v[0], y = v[1], z = v[2];
31
        return x*x*y + std::sin(y) + x*std::exp(z);
32
    };
33
34
    std::vector<double> x = {1.2, -0.7, 0.3};
35
36
    for (std::size_t i = 0; i < x.size(); ++i) {
37
        double dfi = partial_derivative(f, x, i);
38
        std::cout << "Partial derivative wrt x[" << i << "]: " << dfi << "\n";
39
    }
40
41
    return 0;
42
}
43

1
#include <cmath>
2
#include <functional>
3
#include <iostream>
4
#include <limits>
5
#include <vector>
6
7
std::vector<double> gradient(const std::function<double(const std::vector<double>&)>& f,
8
                             const std::vector<double>& x) {
9
    std::size_t n = x.size();
10
    std::vector<double> g(n);
11
    std::vector<double> xt = x;
12
    const double eps = std::numeric_limits<double>::epsilon();
13
14
    for (std::size_t i = 0; i < n; ++i) {
15
        double xi = x[i];
16
        double h = std::sqrt(eps) * std::max(1.0, std::abs(xi));
17
        xt[i] = xi + h; double fp = f(xt);
18
        xt[i] = xi - h; double fm = f(xt);
19
        g[i] = (fp - fm) / (2.0 * h);
20
        xt[i] = xi; // restore
21
    }
22
    return g;
23
}
24
25
int main() {
26
    // Same test function as before for consistency.
27
    auto f = [](const std::vector<double>& v) -> double {
28
        double x = v[0], y = v[1], z = v[2];
29
        return x*x*y + std::sin(y) + x*std::exp(z);
30
    };
31
32
    std::vector<double> x = {1.2, -0.7, 0.3};
33
34
    std::vector<double> g = gradient(f, x);
35
36
    std::cout << "f(x) = " << f(x) << "\n";
37
    std::cout << "Gradient: [" << g[0] << ", " << g[1] << ", " << g[2] << "]\n";
38
39
    return 0;
40
}
41

1
#include <cmath>
2
#include <iostream>
3
#include <vector>
4
5
struct Dual {
6
    double val;                 // primal value
7
    std::vector<double> der;    // derivative vector (size n)
8
9
    Dual() : val(0.0) {}
10
    Dual(double v, std::size_t n) : val(v), der(n, 0.0) {}
11
12
    static Dual variable(double v, std::size_t n, std::size_t idx) {
13
        Dual d(v, n);
14
        d.der[idx] = 1.0;       // seed direction e_idx
15
        return d;
16
    }
17
};
18
19
// Helper: elementwise operations on derivative vectors
20
static inline std::vector<double> add(const std::vector<double>& a, const std::vector<double>& b) {
21
    std::vector<double> r(a.size());
22
    for (std::size_t i = 0; i < a.size(); ++i) r[i] = a[i] + b[i];
23
    return r;
24
}
25
static inline std::vector<double> sub(const std::vector<double>& a, const std::vector<double>& b) {
26
    std::vector<double> r(a.size());
27
    for (std::size_t i = 0; i < a.size(); ++i) r[i] = a[i] - b[i];
28
    return r;
29
}
30
static inline std::vector<double> scal(const std::vector<double>& a, double s) {
31
    std::vector<double> r(a.size());
32
    for (std::size_t i = 0; i < a.size(); ++i) r[i] = a[i] * s;
33
    return r;
34
}
35
36
// Operators for Dual
37
Dual operator+(const Dual& a, const Dual& b) {
38
    Dual r; r.val = a.val + b.val; r.der = add(a.der, b.der); return r;
39
}
40
Dual operator-(const Dual& a, const Dual& b) {
41
    Dual r; r.val = a.val - b.val; r.der = sub(a.der, b.der); return r;
42
}
43
Dual operator*(const Dual& a, const Dual& b) {
44
    // Product rule: d(a*b) = a'*b + a*b'
45
    Dual r; r.val = a.val * b.val;
46
    // r.der = a.der*b.val + b.der*a.val
47
    std::size_t n = a.der.size();
48
    r.der.assign(n, 0.0);
49
    for (std::size_t i = 0; i < n; ++i) r.der[i] = a.der[i] * b.val + b.der[i] * a.val;
50
    return r;
51
}
52
53
// Dual with double interactions (minimal set used below)
54
Dual operator+(const Dual& a, double b) { Dual r; r.val = a.val + b; r.der = a.der; return r; }
55
Dual operator+(double a, const Dual& b) { return b + a; }
56
Dual operator*(const Dual& a, double b) { Dual r; r.val = a.val * b; r.der = scal(a.der, b); return r; }
57
Dual operator*(double a, const Dual& b) { return b * a; }
58
59
// Elementary functions (extend as needed)
60
Dual sin(const Dual& a) {
61
    Dual r; r.val = std::sin(a.val); r.der = scal(a.der, std::cos(a.val)); return r;
62
}
63
Dual exp(const Dual& a) {
64
    double ev = std::exp(a.val);
65
    Dual r; r.val = ev; r.der = scal(a.der, ev); return r;
66
}
67
68
// Example function f(x,y,z) = x^2*y + sin(y) + x*exp(z)
69
Dual f(const std::vector<Dual>& X) {
70
    const Dual& x = X[0];
71
    const Dual& y = X[1];
72
    const Dual& z = X[2];
73
    return (x * x) * y + sin(y) + x * exp(z);
74
}
75
76
int main() {
77
    std::size_t n = 3;                      // number of variables
78
    std::vector<double> point = {1.2, -0.7, 0.3};
79
80
    // Lift to Dual variables with appropriate seeds
81
    std::vector<Dual> X(n);
82
    for (std::size_t i = 0; i < n; ++i) X[i] = Dual::variable(point[i], n, i);
83
84
    Dual y = f(X);                          // one pass computes value and gradient
85
86
    std::cout << "f(x) = " << y.val << "\n";
87
    std::cout << "grad f = [";
88
    for (std::size_t i = 0; i < n; ++i) {
89
        std::cout << y.der[i] << (i + 1 < n ? ", " : "]\n");
90
    }
91
    return 0;
92
}
93