C++ - Logistic Regression Backpropagation with Gradient Descent

Question

I implemented binary logistic regression for a single datapoint trained with the backpropagation algorithm to calculate derivatives for a gradient descent optimizer.

I am primarily looking for feedback on how I approached the functions that return optional derivatives. I used nullptr as a flag to indicate that derivatives are not needed, but was wondering whether the alternatives would be better (i.e., using a boolean flag, using an optional, using a non-const reference, etc.). The derivatives are optional since inference would not need the derivatives. I have left the forward pass functionality out to keep the code a bit smaller.

I am also interested in performance, style, and correct usage of C++ idioms. I am aware that linear algebraic methods would benefit from vectorized methods. For performance concerns, I'm more interested in things like needless copies and redundant computations rather than vectorization opportunities.

Note that I use the naming convention from Andrew Ng's Deep Learning course of the variable name dVAR in a function f to mean the derivative of f with respect to VAR. For example, in double f(double x, double* dx) dx\$ = \frac{df(x)}{dx}\$.

#include <array>
#include <cmath>

template<int sz>
class Vec {
public:
  Vec<sz>() {}
  Vec<sz>(const std::array<double, sz>& storage) : storage_{storage} {}
  Vec<sz>(std::array<double, sz>&& storage) : storage_{move(storage)} {}
  
  Vec<sz> operator-(const Vec<sz>& other) const {
    std::array<double, sz> result;
    for (size_t i = 0; i < storage_.size(); ++i) {
      result[i] = storage_[i] - other.storage_[i];
    }
    return Vec<sz>{result};
  }
  
  Vec<sz> operator*(const double scalar) const {
    std::array<double, sz> result;
    for (size_t i = 0; i < storage_.size(); ++i) {
      result[i] = scalar * storage_[i];
    }
    return Vec<sz>{result};
  }
  
  double operator[](const size_t idx) const {
    return storage_[idx];
  }
  
  bool operator==(const Vec& other) const {
    return storage_ == other.storage_;
  }
  
  size_t size() const {
    return sz;
  }
private:
  std::array<double, sz> storage_;
};

template<int sz>
Vec<sz> operator*(const double scalar, Vec<sz> v) {
  return v * scalar;
}

// Computes the sigmoid function and optionally its derivative.
double sigmoid(const double x, double* dx=nullptr) {
  const double ret = 1 / (1 + std::exp(-x));
  if (dx) {
    *dx = ret * (1 - ret);
  }
  return ret;
}

// Compute the cross entropy and optionally its derivative with respect to y_hat.
double cross_entropy(const int y, const double& y_hat, double* dy_hat=nullptr) {
  const double ret = - (y * std::log(y_hat) + (1 - y) * std::log(1 - y_hat));
  
  if (dy_hat) {
    *dy_hat = (y / (1 - y_hat)) + ((1 - y) / y_hat);
  }
  
  return ret;
};

// Compute the dot product and optionally its derivative with respect to w.
template<int size>
double dot_product(const Vec<size>& w,
                   const Vec<size>& x,
                   Vec<size>* dw=nullptr) {
  double ret = 0.0;
  for (size_t i = 0; i < w.size(); ++i) {
    ret += w[i] * x[i];
  }
  
  if (dw) {
    *dw = x;
  }
  return ret;
}

// Do gradient descent of logistic regression updating weights and bias in place.
template<int size>
double logistic_regression(Vec<size>& w,
                           double& b,
                           const Vec<size>& x,
                           const int y,
                           const int num_epochs=100,
                           const double learning_rate=0.001) {
  double final_loss;
  for (int i = 0; i < num_epochs; ++i) {
    // Derivatives.
    Vec<size> dw;
    double dz;
    double da;
    
    // Forward pass.
    const double z = dot_product(w, x, &dw) + b;
    const double a = sigmoid(z, &dz);
    const double loss = cross_entropy(y, a, &da);
    
    // Backward pass.
    w = w - learning_rate * (dw * dz * da);
    b = b - learning_rate * (dz * da);
    
    final_loss = loss;
  }
  return final_loss;
}

Answer 1

About the Vec class: Every function in this class could be constexpr. (Especially size().) At least operator== and size() could also be noexcept. constexpr and noexcept could vastly improve code generation in most cases.

(Not all of the functions should really be noexcept, though. operator[] probably shouldn’t be, for example, because it might fail for out-of-range indexes. That doesn’t mean you should literally throw an exception! A non-noexcept function doesn’t have to throw. You’re just leaving it as an option, maybe just for in debug mode (and in release mode, remove all exceptions for max performance). You should read noexcept as “this function can’t possibly fail” and anything else as “this function can fail, and might throw when it does (or it might not)”.)

Vec<sz>(const std::array<double, sz>& storage) : storage_{storage} {}
Vec<sz>(std::array<double, sz>&& storage) : storage_{move(storage)} {}

There’s really no point in the second constructor. Arrays are not moveable. The contents of an array might be moveable, but in this case the contents are doubles, and those are also not moveable. So you have a copy-only container of copy-only stuff… the moving constructor above ends up doing exactly the same thing as the copying constructor, making it pointless.

Vec<sz> operator-(const Vec<sz>& other) const {

It’s generally a good idea to make your binary operations non-member functions—it’s generally a good idea to make as many things as possible non-member functions, for the sake of encapsulation.

The usual way to do this is to define the assignment versions of binary ops in-class, and then write the regular binary ops in terms of the assignment versions. Or in plain C++lish:

class Type
{
public:
    auto operator+=(Type const& t) -> Type&
    {
        // define += operation...
        return *this;
    }

    // rest of class...
};

auto operator+(Type a, Type const& b) -> Type
{
    return a += b;
}

So you should probably write an operator-= for Vec<sz>, and then define operator- outside of the class in terms of operator-=. (You probably want an operator-= anyway, for reasons I’ll explain later.)

Vec<sz> operator-(const Vec<sz>& other) const {
  std::array<double, sz> result;
  for (size_t i = 0; i < storage_.size(); ++i) {
    result[i] = storage_[i] - other.storage_[i];
  }
  return Vec<sz>{result};
}

result seems to be completely superfluous here. You store the calculation results in result, and then copy them into the Vec<sz>’s internal storage_… why not skip that step, and store the calculation results directly into the storage_? Like so:

Vec<sz> operator-(const Vec<sz>& other) const {
  auto result = Vec<sz>{};
  for (size_t i = 0; i < storage_.size(); ++i) {
    result.storage_[i] = storage_[i] - other.storage_[i];
  }
  return result;
}

With C++17’s guaranteed return value elision, you can end up with no copying at all.

Vec<sz> operator*(const double scalar) const {

This is another function that would be better as a non-member. In this case, you also have symmetry going on:

template <int sz>
auto operator*(Vec<sz> const& v, double scalar) -> Vec<sz> { ... }
template <int sz>
auto operator*(double scalar, Vec<sz> const& v) -> Vec<sz> { ... }

You could make one or both of those a friend of the class, or define both in terms of a member function (perhaps named scale()).

operator== should also be a non-member, and it should come along with operator!=. (operator[] has to be a member, of course.)

double sigmoid(const double x, double* dx=nullptr) {

I’m not a fan of default arguments. They make writing functions easier, but using them harder… which is the worst kind of feature, because you only write a function once, but you use them many times.

Even worse, in this case you’re actually ruining the function by making it less efficient for every possibly usage. If I don’t want the the derivative, I still have to pay for the extra function parameter and the check. If I do, I’m paying for a mostly unnecessary check (branching is slow!).

And on top of all that, if you don’t try to cram all possible usage scenarios into a single function, you avoid the entire question of whether to use nullptr or std::optional, and make the function easier to use because you can return the derivative directly:

double sigmoid(const double x) {
  return 1 / (1 + std::exp(-x));
}

std::tuple<double, double> sigmoid_with_derivative(const double x) {
  const double ret = 1 / (1 + std::exp(-x));
  return {ret, ret * (1 - ret)};
}

// usage:
auto a = sigmoid(x); // when I don't want the derivative
auto [a, dx] = sigmoid_with_derivative(x); // when I do want it

The same idea applies to cross_entropy() and dot_product().

double cross_entropy(const int y, const double& y_hat, double* dy_hat=nullptr) {

I’m guessing the & here is a typo, but just in case not: there’s (usually!) nothing to be gained by passing built-in types like double by reference.

template<int size>
double logistic_regression(Vec<size>& w,
                           double& b,
                           const Vec<size>& x,
                           const int y,
                           const int num_epochs=100,
                           const double learning_rate=0.001) {

Hm… okay, the last two arguments here are probably legitimate use case for default arguments.

As for the first two arguments… I’m really not a fan of “out” arguments, if that’s what w and b are intended to be. But I don’t know what the intended usage of this function looks like, so I don’t know if maybe this is one of the very, very rare cases where they make sense.

Assuming that isn’t the case, a better form for this function might be something like:

template<int size>
std::tuple<double, Vec<size>, double>
logistic_regression(Vec<size> w,
                           double b,
                           const Vec<size>& x,
                           const int y,
                           const int num_epochs=100,
                           const double learning_rate=0.001) {

  // ... everything else in the function is the same except the last line...

  return {final_loss, w, b};
}

or:

template<int size>
std::tuple<double, Vec<size>, double>
logistic_regression(Vec<size> const w_,
                           double const b_,
                           const Vec<size>& x,
                           const int y,
                           const int num_epochs=100,
                           const double learning_rate=0.001) {
  auto ret = std::tuple{0.0, w_, b_};
  auto&& final_loss = std::get<0>(ret);
  auto&& w = std::get<1>(ret);
  auto&& b = std::get<2>(ret);

  // ... everything else in the function is the same except the last line...

  return ret;
}

Just a few final comments:

It’s good that you’re not fussing over vectorization, because, really, compilers are smart enough to already vectorize your loops. Even in the case that you’re working with a compiler that’s not, it’s pretty trivial to signal you want a loop vectorized with intrinsics or OMP or other means.

You are correct worry that your biggest performance drain will probably be poor algorithms and/or unnecessary copying. I don’t see any obvious signs of redundant or unnecessary calculations, so you might be okay on that front. That leaves unnecessary copying.

You’re using arrays for your data types, rather than, say, vectors. That can actually be a performance pessimization, especially if sz is large. Arrays are non-moveable; they can only be copied. Vectors can be copied or moved, and moves are so fast they’re basically free. On the other hand, when vectors do need to be copied, the cost can be many times more than the cost of copying an array, due to the overhead of allocation.

Note that I’m not suggesting you should use vectors instead! I’m just saying that when using arrays, you can’t count on moves. Every move is actually a copy. So if you are going to use arrays—which is a perfectly sensible thing to do for reasonably small sz—you need to be very careful to avoid both copies and moves.

You note (via a tag) that you’re using C++17. That’s good! That works out a lot in your favour, because C++17 brought in guaranteed elision. What does that mean?

It means that for a function like this:

template <int size>
auto dot_product_with_derivative(
        Vec<size> const& w,
        Vec<size> const& x)
    -> std::tuple<double, Vec<size>>
{
    auto ret = std::tuple{0.0, x};

    for (std::size_t i = 0; i < size; ++i)
        ret += w[i] * x[i];

    return ret;
}

used like this:

auto [z, dw] = dot_product_with_derivative(w, x);

there is only a single array copy done (when the array in x is copied into what will eventually be dw). All other temporary/intermediate values are elided away. This is guaranteed by the standard starting C++17.

The only place I can see where you (might!) have a superfluous copy is in the line:

w = w - learning_rate * (dw * dz * da);

That’s because (dw * dz * da) has to create a temporary Vec<size> (let’s call it tmp1)…

… and then learning_rate * tmp1 has to create another temporary Vec<size> (let’s call it tmp2)…

… and then w - tmp2 has to create ANOTHER temporary…

… when then gets assigned into w.

That’s THREE temporary Vec<size> objects, each of which requires copying an array of doubles.

A smart compiler might be able to elide away that third temporary… but why risk it? Why not write operator-= and do:

w -= learning_rate * (dw * dz * da);

But that still leaves 2 temporaries.

Is there any reason you want to force (dw * dz * da) to be calculated first? That creates a temporary Vec<size> that you then multiply with learning_rate.

Why not reorder your calculation as: (learning_rate * dz * da) * dz?

w -= (learning_rate * dz * da) * dz;

Now (learning_rate * dz * da) only creates a temporary double… not a temporary Vec<size>. A temporary double is basically free.

Then (learning_rate * dz * da) * dz creates a temporary Vec<size>. That’s unavoidable. Let’s call that tmp1.

And then w -= tmp1 creates no temporaries.

So you’ve gone from:

// Up to 3 temporaries created:
w = w - learning_rate * (dw * dz * da);

to:

// Only 1 temporary created at maximum, guaranteed:
w -= (learning_rate * dz * da) * dz;

without changing the actual calculation.

Now, I’m not sure if a smart compiler could do that optimization on its own. I’m not sure how blasé optimizers these days are about ignoring your parentheses, or reordering built-in operations around user-defined operations. Your original line of code might be compiled without all those temporaries. (Probably will, to be honest.) But why risk it? Rewriting x = x - y as x -= y is an optimization as old as time—it’s why operator-= exists in the first place. And reordering your computation so all the double calculations get done first and then only a single Vec<size> calculation has to be done is also a trivial thing to do. Even if it ultimately makes no difference (because optimizing compilers are so smart these days), it doesn’t really cost anything.

Without knowing the expected values for the sz parameter of Vec<sz>, I can’t comment on whether using arrays is a good idea. If sz is fairly small, then sure. Even if sz ranges into a hundred or so, maybe, if you’re not making a ton of these objects. But past a certain point, trucking the arrays around will get so unwieldy that it might make more sense to use vectors. With vectors, you would pay a lot more at construction time, and if copying… but you get super-fast moves. Might be worth it, especially if the rest of the code is written to take advantage of moving.