# Understanding the stochastic average gradient (SAG) algorithm used in sklearn

For pedagogical purposes I've been trying to create my own implementation of the stochastic average gradient (SAG) algorithm in a logistic regression framework. Page 10 of the associated paper describes the algorithm for SAG, which seems relatively straightforward.

$$y_{i}$$ stores the gradient of the loss function with respect to the weights for randomly selected observation $$i \in \{1, ..., n\}$$ (e.g. if there are p features, the gradient will be $$p$$ x 1 or 1 x $$p$$). $$y$$ thus is of dimension $$n$$ x $$p$$ or $$p$$ x $$n$$.

$$d$$ is the sum of the gradients, $$\sum_{i=1}^{n}y_{i}$$

$$k$$ is the number of iterations.

$$\eta$$ is the learning rate

Thus the algorithm with respect to a single weight $$x$$ is defined as:

d = 0
y_i = 0 for i = 1, ..., n

for k = 0, ..., max_iterations
Sample i from {1, ..., n}
d = d - y_i + f'_i(x)
y_i = f'_i(x)
x = x - (eta / n) * d
end


In other words, we initialize a table of loss gradients with respect to the weights for each observation at $$0$$. Likewise, the sum of these weights across all observations will start at $$0$$. Then we randomly sample one observation, compute the loss gradient, and update the sum of gradients by adding this gradient and subtracting the previous gradient for that same sampled observation. We then subtract the product of the sum of the gradients (divided by the total number of observations) and the learning rate from our existing weights to obtain the new weights.

It feels like I'm missing something here, since by initializing at 0, we're effectively slowly building a stochastic average that is aggressively weighed down by the fact that we average across all observations (divide by $$n$$) even before we've generated a non-zero gradient. The first iteration would in effect be equal to stochastic GD, except dividing by $$n$$ observations, for example. My own implementation in Python vastly underperforms sklearn, and standard stochastic gradient descent for that matter.

Am I misinterpreting how the gradients are generated or updated?

1. Instead of initializing y_i =0, instead spend one pass over the data and initialize y_i = f'_i(x_0)
2. The more practical fix is the do one epoch SGD over the shuffled data, and record the gradient y_i = f'_i(x_i). After the first epoch, then switch to SAG or SAGA.