|
| 1 | +/* |
| 2 | +This demo implements logistic regression with CUDA Thrust and clad automatic |
| 3 | +differentiation. It differentiates a single-example logistic loss, then runs SGD |
| 4 | +on a small four-document batch by invoking an L2-regularized two-document loss |
| 5 | +twice. The program prints the single-example loss and gradient, reports |
| 6 | +training loss periodically, and finally outputs the batch accuracy. |
| 7 | +*/ |
| 8 | + |
| 9 | +#include "clad/Differentiator/Differentiator.h" |
| 10 | +#include "clad/Differentiator/ThrustBuiltins.h" |
| 11 | +#include "clad/Differentiator/ThrustDerivatives.h" |
| 12 | + |
| 13 | +#include <thrust/device_vector.h> |
| 14 | +#include <thrust/fill.h> |
| 15 | +#include <thrust/functional.h> |
| 16 | +#include <thrust/host_vector.h> |
| 17 | +#include <thrust/inner_product.h> |
| 18 | +#include <thrust/iterator/permutation_iterator.h> |
| 19 | +#include <thrust/reduce.h> |
| 20 | +#include <thrust/transform.h> |
| 21 | + |
| 22 | +#include <cmath> |
| 23 | +#include <iostream> |
| 24 | + |
| 25 | +using Vec = thrust::device_vector<double>; |
| 26 | +static inline double dot(const Vec& a, const Vec& b) { |
| 27 | + return thrust::inner_product(a.begin(), a.end(), b.begin(), 0.0); |
| 28 | +} |
| 29 | +static inline double sigmoid(double z) { return 1.0 / (1.0 + std::exp(-z)); } |
| 30 | + |
| 31 | +// Minimal single-document logistic loss. |
| 32 | +// x: (V), w: (V), y in {0,1} |
| 33 | +double logistic_loss_single(const Vec& x, const Vec& w, double b, double y) { |
| 34 | + double logit = dot(x, w) + b; |
| 35 | + double p = sigmoid(logit); |
| 36 | + const double eps = 1e-9; |
| 37 | + return -y * std::log(p + eps) - (1.0 - y) * std::log(1.0 - p + eps); |
| 38 | +} |
| 39 | + |
| 40 | +// L2-regularized 2-doc batch loss (0.5*(l0+l1) + 0.5*lambda*(||w||^2 + b^2)) |
| 41 | +double logistic_loss_batch2_prepared_l2(const Vec& x0, const Vec& x1, |
| 42 | + const Vec& w, double b, double y0, |
| 43 | + double y1, double lambda) { |
| 44 | + double l0 = logistic_loss_single(x0, w, b, y0); |
| 45 | + double l1 = logistic_loss_single(x1, w, b, y1); |
| 46 | + double w2 = dot(w, w); |
| 47 | + return 0.5 * (l0 + l1) + 0.5 * lambda * w2 + 0.5 * lambda * (b * b); |
| 48 | +} |
| 49 | + |
| 50 | +static inline void zero(Vec& v) { thrust::fill(v.begin(), v.end(), 0.0); } |
| 51 | + |
| 52 | +static inline double predict_prob(const Vec& xi, const Vec& w, double b) { |
| 53 | + return sigmoid(dot(xi, w) + b); |
| 54 | +} |
| 55 | + |
| 56 | +static inline void sgd_step(Vec& w, const Vec& dw, double& b, double db, |
| 57 | + double lr) { |
| 58 | + thrust::transform( |
| 59 | + w.begin(), w.end(), dw.begin(), w.begin(), |
| 60 | + [=] __device__(double wi, double gi) { return wi - lr * gi; }); |
| 61 | + b -= lr * db; |
| 62 | +} |
| 63 | + |
| 64 | +int main() { |
| 65 | + std::cout << "Running minimal logistic regression demo." << std::endl; |
| 66 | + |
| 67 | + const int V = 8; |
| 68 | + // Dense feature vector and weights. |
| 69 | + double hx[] = {2, 0, 1, 0, 0, 0, 1, 0}; |
| 70 | + Vec x(hx, hx + V); |
| 71 | + Vec w(V, 0.1); |
| 72 | + double b = 0.0; |
| 73 | + double y = 1.0; |
| 74 | + |
| 75 | + auto grad = clad::gradient(logistic_loss_single); |
| 76 | + Vec dx(V), dw(V); |
| 77 | + double db = 0.0; |
| 78 | + double dy = 0.0; |
| 79 | + |
| 80 | + grad.execute(x, w, b, y, &dx, &dw, &db, &dy); |
| 81 | + double loss = logistic_loss_single(x, w, b, y); |
| 82 | + |
| 83 | + thrust::host_vector<double> hdw = dw; |
| 84 | + std::cout << "Loss: " << loss << "\nGradient wrt w: "; |
| 85 | + for (int i = 0; i < V; ++i) |
| 86 | + std::cout << hdw[i] << " "; |
| 87 | + std::cout << std::endl; |
| 88 | + |
| 89 | + // 4-document batch demo. |
| 90 | + { |
| 91 | + std::cout << "\nRunning SGD on 4-doc batch..." << std::endl; |
| 92 | + double hX4[] = { |
| 93 | + 2, 0, 1, 0, 0, 0, 1, 0, // doc 0 (y=1) |
| 94 | + 0, 3, 0, 1, 0, 0, 0, 0, // doc 1 (y=0) |
| 95 | + 0, 1, 2, 0, 0, 0, 0, 0, // doc 2 (y=1) |
| 96 | + 0, 0, 0, 0, 2, 1, 0, 0 // doc 3 (y=0) |
| 97 | + }; |
| 98 | + Vec X4(hX4, hX4 + 4 * V); |
| 99 | + Vec a0(X4.begin() + 0 * V, X4.begin() + 1 * V); |
| 100 | + Vec a1(X4.begin() + 1 * V, X4.begin() + 2 * V); |
| 101 | + Vec a2(X4.begin() + 2 * V, X4.begin() + 3 * V); |
| 102 | + Vec a3(X4.begin() + 3 * V, X4.begin() + 4 * V); |
| 103 | + double y_a0 = 1.0, y_a1 = 0.0, y_a2 = 1.0, y_a3 = 0.0; |
| 104 | + |
| 105 | + auto loss2 = clad::gradient(logistic_loss_batch2_prepared_l2); |
| 106 | + Vec dA0(V), dA1(V), dA2(V), dA3(V), dW4(V); |
| 107 | + double db4 = 0.0, dlam4 = 0.0; |
| 108 | + const int iters4 = 50; |
| 109 | + const double lr4 = 0.1; |
| 110 | + const double lambda4 = 1e-2; |
| 111 | + for (int t = 0; t < iters4; ++t) { |
| 112 | + zero(dW4); |
| 113 | + zero(dA0); |
| 114 | + zero(dA1); |
| 115 | + zero(dA2); |
| 116 | + zero(dA3); |
| 117 | + db4 = dlam4 = 0.0; |
| 118 | + |
| 119 | + double dy_dummy = 0.0; |
| 120 | + loss2.execute(a0, a1, w, b, y_a0, y_a1, lambda4, &dA0, &dA1, &dW4, &db4, |
| 121 | + &dy_dummy, &dy_dummy, &dlam4); |
| 122 | + loss2.execute(a2, a3, w, b, y_a2, y_a3, lambda4, &dA2, &dA3, &dW4, &db4, |
| 123 | + &dy_dummy, &dy_dummy, &dlam4); |
| 124 | + |
| 125 | + thrust::transform(dW4.begin(), dW4.end(), dW4.begin(), |
| 126 | + [] __device__(double g) { return 0.5 * g; }); |
| 127 | + db4 *= 0.5; |
| 128 | + |
| 129 | + sgd_step(w, dW4, b, db4, lr4); |
| 130 | + |
| 131 | + if ((t % 10) == 0 || t == iters4 - 1) { |
| 132 | + double Lpair1 = |
| 133 | + logistic_loss_batch2_prepared_l2(a0, a1, w, b, y_a0, y_a1, lambda4); |
| 134 | + double Lpair2 = |
| 135 | + logistic_loss_batch2_prepared_l2(a2, a3, w, b, y_a2, y_a3, lambda4); |
| 136 | + double L4 = 0.5 * (Lpair1 + Lpair2); |
| 137 | + std::cout << "iter4 " << t << ": loss4=" << L4 << std::endl; |
| 138 | + } |
| 139 | + } |
| 140 | + |
| 141 | + int ok = 0; |
| 142 | + ok += ((predict_prob(a0, w, b) >= 0.5) == (y_a0 == 1.0)); |
| 143 | + ok += ((predict_prob(a1, w, b) >= 0.5) == (y_a1 == 1.0)); |
| 144 | + ok += ((predict_prob(a2, w, b) >= 0.5) == (y_a2 == 1.0)); |
| 145 | + ok += ((predict_prob(a3, w, b) >= 0.5) == (y_a3 == 1.0)); |
| 146 | + std::cout << "Batch-4 accuracy: " << (double)ok / 4.0 << std::endl; |
| 147 | + } |
| 148 | + |
| 149 | + return 0; |
| 150 | +} |
0 commit comments