Machine Learning 1. Lecture 2. Stochastic Gradient and the Linear Classifier

What this lecture is about

A dense lecture. Today: linear models of classification and regression and the stochastic gradient method; with these I will lead you up to the problem statements and methods used in training neural networks, and we will reach the networks themselves in the next lecture. At the end there will be one more important topic — the probabilistic statement of the classification problem and the probabilistic meaning of regularization.

Four steps: regression, binary and multiclass classification

In the previous lecture I showed you the method of least squares, and we looked at how linear, and in particular polynomial, regression works. Let us once more bring all this into a scheme of four steps and then extend it to classification.

The regression problem. We are given a training sample — «object — answer» pairs $(x_i,y_i{)}_{i=1}^{\ell }$, the answers are real numbers, $y_i\in \mathbb{R}$. We solve it like this:

1. Fix the model — here the linear regression model, a linear combination of all features with coefficients $w$:
   $$a(x,w)=⟨w,x⟩=\sum _{j=1}^n{w}_j{f}_j(x).$$
   Get used to writing it as a scalar product: from now on I will mostly use this vector notation.
2. Fix the loss function — let it be quadratic, so that we get the method of least squares.
3. Reduce learning to optimization: sum the losses over all the objects of the training sample and fit the parameter vector $w$ so that on the given objects the model gives answers as close to the correct ones as possible:
   $$Q(w)=\sum _{i=1}^{\ell }(⟨w,x_i⟩-y_i{)}^2\to \underset{w}{\min }.$$
4. Check on test data. For machine learning, optimization is only a tool. We must build a model that has predictive power, or generalization ability — people say it differently. We check this on a test sample, on data not used in training at all: the mean loss on it must also be adequate, small enough.

The binary classification problem. Now the same scheme for classification; I deliberately made two very similar slides. The answers are classes, and there are two classes: $y_i\in \{-1,+1\}$. The problem is called binary because we make a discrete «yes» or «no» decision: treat or not treat, grant a loan or not, is there oil in a region or not. As we saw in the previous lecture, there are very many such applied problems.

The model is again a scalar product, but now we are interested not in its value, only in its sign:

If the scalar product is greater than zero, the object $x$ goes to class $+1$; if less, to class $-1$.

How should the loss function be built? It would seem nothing could be simpler: take the condition that the classifier errs. The classifier returns $\pm 1$, the correct answer $y$ is also $\pm 1$. If their product is negative, they have different signs, that is, the answer $a$ the classifier gives does not coincide with the correct answer $y$. The function $\mathrm{sign}$ does not change the sign of its argument, so the scalar product itself can be substituted for $a$, and an amusing and very useful construction arises:

What is its geometric meaning? The equation $⟨w,x⟩=0$ defines a hyperplane in ${\mathbb{R}}^n$ with direction vector $w$, which divides the space into two half-spaces. In the positive direction, where $w$ points, lie the objects of class $+1$; in the opposite half, the objects of class $-1$. The quantity $⟨w,x⟩$ shows, roughly speaking, the signed distance from the object $x$ to the separating hyperplane. If $x$ is an object of class $+1$ and the angle between it and $w$ is acute, it has landed in the right half of the space, and the product $⟨w,x⟩y$ is greater than zero. If it is less than zero, the object has landed in the wrong half, and an error has occurred on it.

The threshold loss — these square brackets — is unpleasant because its values are discrete. If we build the functional as a sum of discrete losses, we get a discrete-valued criterion, very inconvenient for optimization. So we will minimize not it but an upper bound on it: take a continuous, or better a smooth, function $L(M)$ that bounds the threshold one from above:

Then we approximately minimize the number of errors on the training sample, and we will be able to get a fairly decent vector $w$ this way. Here is a very important principle, which can now be called the dominant one in classification theory — the principle of approximating the threshold loss function. We have not only made the problem nicer for numerical methods: in this way we can also penalize an object for approaching the separating hyperplane; more on that a bit later. The fourth step, checking on test data, we of course do not forget.

The multiclass classification problem. What if there are many classes? Now $Y$ is a finite set of classes, in general more than two. We would still like to keep the model linear. This construction helps: again we compute a scalar product, but each class has its own direction vector, and we assign the object to the class for which the scalar product is largest:

How is this to be interpreted? Each vector says: the objects of my class lie in that direction, and the farther you go there, the more objects of my class there are. Since each class has its own preferred direction, somewhere at the junction the models have to sort things out among themselves and come to an agreement; most likely there will be errors there, and we need a criterion that makes the errors as few as possible.

The loss function is again binary: if we guessed the class wrong, we must be penalized. But an object $x$ of class $y$ can be confused with different classes, and confusing it with one class is not as bad as with ten. So it is proposed to sum the errors over all the wrong classes $z\ne y$: the scalar product $⟨w_y,x⟩$ should have been the largest among all classes, and the number of classes for which it turned out smaller is the loss. The binary loss is still inconvenient, but we already know how to deal with it: turn the logical condition into a margin, now for each pair «correct class $y$ — wrong class $z$»:

The larger it is, the better; if the scalar product for the correct class is larger than for all the others, all such margins are positive and there is no error on the object. The learning method is still empirical risk minimization, only one more sum appears: over all the objects of the sample and, for each object, over all the wrong classes:

Here we use the «one against one» principle: all classes are compared with each other pairwise, and each pair of classes is handled as a binary classification. There are other strategies for multiclass classification too.

The separating surface and the margin

Does it matter so much for this scheme that the model is linear? Of course not: that was just an illustration to make things easier for you. Nonlinear separating surfaces in classification? Easy: the scalar product will have to be replaced by some other function $g(x,w)$ of the object and the parameter vector. And the vector $w$ is then not at all obliged to have the same dimension as $x$ — it all depends on how the model is built.

If we have written out the discriminant function and found the parameter vector $w$, we know what the separating surface between the classes looks like: it is the locus of points $x$ of the object space where $g(x,w)=0$. The margin is still a kind of signed distance to the separating surface, but now the surface can be curved, anything at all, depending on how the function $g$ is defined. The meaning is the same: the margin is negative if and only if the algorithm makes an error on the object.

Error analysis by margins

Most classification methods rest on introducing something like a discriminant function, and here is a very useful plot that can tell a lot about how well the problem has been solved. On the training sample, and on the test one too — as long as the correct answer is known — we can compute the margins of all objects, sort them in increasing order and look at the result as a plot.

On the right are the objects classified very reliably: they lie deep inside their class, all is well with them. Where the margin is negative there are errors, but errors differ. Errors along the boundary between classes may happen under the influence of random factors. But objects with a large negative margin most likely fit our classification model very badly. This may be due to noise: for example, some features were measured completely wrong, or even the answer $y$ was mixed up. There can be very many reasons why an object flies off into the noise region, and it is always useful to look at this plot. It is even more useful to overlay the same plot for the test sample: if on the test the curve shifts too far to the right, that is, the margins there are much worse, that is the effect of overfitting.

The multiclass case

In multiclass classification a separate discriminant function has to be built for each class:

The separating surface can still be defined, but now between each pair of classes: it is the locus of points where two discriminant functions, $g_y$ and $g_z$, are equal. The margin is first defined relative to each wrong class,

and taking the minimum over all wrong classes we get what is called the margin of the object:

It is negative if and only if the algorithm errs on the object. The criterion is the minimum of the approximated empirical risk:

In exactly the same way, the margins in the multiclass case can be sorted to see how many objects are noise and how many are borderline, where the classification is uncertain: such a plot gives a great deal of information for error analysis.

Continuous approximations of the threshold loss

A question I have so far put off and swept under the carpet: what continuous approximations of the threshold loss are there? There are a heap of spoilers here: we will study all these methods in the course, and I will explain all the names.

You will say that some of these functions are not upper approximations at all: how is $H$ an upper approximation? But we can add one to it, and it will be an upper approximation, and that will not affect the solution of the optimization problem. And why a quadratic function that also penalizes overly large positive margins — we will get to that miracle in due course, towards the end of the semester, and see when it happens.

The linear classifier is a model of a neuron

Now let us talk in more detail about binary classification, and specifically about the linear classifier model. It is a model of a neuron, the elementary building block of all modern neural networks, but also a perfectly useful classification model in its own right, which works in many applications. It was invented long ago, in 1943.

A nerve cell is of course much, much more complex than any picture, but this is a model, and a very physical one. A nerve cell has many input processes — dendrites, which can number in the thousands. They let charge into the cell; once inside in a significant amount, exceeding some activation threshold, the charge activates the cell, it becomes excited and transmits an impulse along a long axon to the next cell, where the axon is connected to the dendrites of one or more other cells. So the impulse is passed from cell to cell. What matters for us in this model? Charge, as we know from physics, is an additive quantity: charges add up. Dendrites may let through not all the charge but part of it, and may even inhibit its passage. Hence the linear model: $f_j$ are the charges at the input dendrites; with weights $w_j$ they get inside, form a total charge, which is compared with the activation threshold $w_0$, and if activation occurs, the function $\sigma$ — for example, a threshold one — emits an impulse.

Linear models, first, are what multilayer neural networks are built from, and second, the linear model as a classifier on its own is itself used to solve applied problems.

The stochastic gradient method

We agreed that we use some loss function, and for generality I will write it as ${\mathcal{L}}_i(w)=\mathcal{L}(w,x_i)$. In regression the loss depends on the object $x_i$, the correct answer $y_i$ and the parameter vector $w$; in classification, the same. So as not to be distracted by secondary details of whether we are solving regression or classification, let us assume that for each object some loss function is given, smooth — that is important — and depending on the parameters $w$, and we need to minimize the sum of such functions over all the objects of the sample. In each particular problem ${\mathcal{L}}_i$ can be written as a formula; the main thing is that the formula will have to be differentiated.

Gradient descent

For some this is a reminder, for others fresh news — how numerical optimization methods work. I will not go deep into details: we need only the simplest method, the first-order gradient method. Any schoolchild can understand it. Take some point as the initial approximation, and then go in the direction of the fastest decrease of the function. How do we know which direction that is? The gradient — the vector of partial derivatives of the functional being optimized — shows the direction of fastest increase of the function, and we need the steepest decrease, so the antigradient:

Each iteration is a little step of size $h$ in the direction of the antigradient. In optimization $h$ is called the gradient step, in machine learning the learning rate. If we substitute our functional here, we get a sum of the gradients of the losses over all the objects of the sample:

This is a very important construction for us, and it has a drawback that immediately leads to an idea. We compute the gradient of the loss at the point $w$ in a loop over all the objects of the sample: for the first object, for the second, and so on for all, and then add them up and get the next iteration. The question: since we have already computed the gradient step for one term, for the $i$th object, can we not shift at once and update the vector $w$? Will convergence then not be faster? Another question arises immediately: must we update one object at a time — perhaps collect a whole pack of objects, a batch, and take a step on it? That is possible too, and it is done.

The idea of speeding up convergence is rather old — it is the Robbins–Monro method, invented back in the fifties. If the original criterion is a sum of a large number of terms, then, since taking a gradient is an additive operation, let us take small gradient steps on each term and update the parameter vector more often.

The algorithm

How to initialize the weights I will put off for now. We need some estimate of the functional already at the zeroth iteration, at least a rough approximation, and there is no point computing it over the whole sample: we are doing everything roughly now, so we take a random subset.

An important point is that the object at each iteration is chosen at random, hence the name: the stochastic gradient method. Why randomness here? It was found that if you keep going through the objects of the sample in a loop in the same order, sometimes the process starts cycling. If the objects are taken in random order, all these artefacts are destroyed and there is no cycling. This is an important discovery, and it is in fact the merit of Robbins and Monro.

We repeat the steps until the process converges. How will we know? The value of $Q$ has stopped changing significantly, or the weights $w$ have stopped changing significantly. In general, in optimization theory the stopping condition is always a heuristic: you can simply set the number of iterations, which is often done in practice, or you can construct very clever combinations of conditions that stop the process earlier, when continuing it no longer makes sense.

Recursive estimation of the functional

Where did the convex combination formula in step 3.4 come from? We want to track how the functional being optimized changes. We saved a great deal by taking only one term, one object, for the gradient step. But the value of $Q$ is computed over the whole sample, and if we recomputed it over the whole sample each time to track $Q$, that would kill all the performance just gained by the stochastic gradient method. We can save by using a recursive formula.

Let us start with the arithmetic mean of the values ${\epsilon }_m,{\epsilon }_{m-1},{\epsilon }_{m-2},\dots$ — at the current step, the previous one, the one before that, and so on. Recursively it is computed very simply — as a convex combination of the previous mean and the last term:

Now replace $\frac{1}{m}$ with $\lambda$. The number of iterations $m$ grows at each step, but we use a trick: we fix $\lambda$. The recursive formula for the arithmetic mean turns into a recursive formula for the exponential moving average:

The coefficients decrease into the past as a geometric progression: information about past errors is forgotten at the rate of a geometric progression. The parameter $\lambda$ is the forgetting rate of the history of the series. A very useful technique; we will come back to it more than once. Comparing the two recursive formulas shows how to set $\lambda$: it is of the order of $\frac{1}{m}$, where $m$ is the number of values we want to average over. If we want to average over roughly the last ten iterations, $\lambda =0.1$; over a hundred, $\lambda =0.01$. In reality the averaging is over all iterations from the very beginning, but roughly speaking one can say it is mainly over the last $\frac{1}{\lambda }$.

Speeding up convergence: momentum and the accelerated gradient

We will talk about this in more detail in the lecture on neural networks; for now, briefly. If we head for the minimum from some random point, at each stochastic gradient step the object is chosen at random and only one term is optimized. On average we move in the direction in which the total function $Q(w)$ decreases, but the motion resembles that of a Brownian particle: a random factor keeps pushing us in different directions — bang one way, bang another — though on average we go where we need to.

Can we go where we need to faster? Yes, by averaging the directions each stochastic gradient step gives. The idea is exactly the same as in the exponential moving average — I told you about it just now for a reason; but there we averaged the value of $Q$ in order to estimate it, and here the motive is quite different.

If $\gamma$ is close to one, we average the gradients over many steps and get a direction that shows more accurately where the criterion decreases. The method was invented by Boris Teodorovich Polyak; it is already considered a classic and a very effective way to speed up the convergence of gradient methods substantially.

The word momentum should be translated into Russian not as «момент» but as «импульс» or «инерция»: it is about the inertia with which the process — a ball — rolls into a potential well, in the direction in which the functional decreases. The oldest gradients we remember with smaller weights, but in fact all gradients from the very first iteration are averaged. Each object adds its gradient vector to the momentum, and this contribution keeps its influence in the vector $v$ for a long time. So we move by inertia in the direction set by all the previous objects.

Another interesting way of speeding things up was invented by Nesterov in 1983 — Nesterov's accelerated gradient. We know the direction in which we are moving by inertia, yet for some reason we compute the gradient at the point where we were. Let us compute it at the point we will approximately get to. It turns out that such a small modification — the gradient's argument changed just a little — also leads to acceleration. Together these two techniques can help the iterative process substantially.

Heuristics of the stochastic gradient method

Let us gallop through various heuristics invented to speed up the stochastic gradient method. Some of them are tied closely to linear models, some are more universal.

Weight initialization

How should the weights be initialized at the first step so as to land at a reasonably good point?

• Small random values. Forget this recommendation as a universal one: it comes from old textbooks, nobody does that any more. It is needed only when the loss functions or the activation functions in neural networks have horizontal asymptotes. Then large weight values may knock us into a region where the derivative of the loss is close to zero and the gradient steps come out small. For such problems — yes, random initialization near zero; but that is far from needed everywhere.
• By correlation. An interesting recommendation for linear models:
  $$w_j:=\frac{⟨y,f_j⟩}{⟨f_j,f_j⟩},f_j=(f_j(x_i){)}_{i=1}^{\ell },y=(y_i{)}_{i=1}^{\ell }.$$
  This is something like the correlation of the vector of values of the $j$th feature on all the objects of the sample with the vector of answers: if a feature is a good predictor for a one-dimensional prediction of $y$, it gets a large weight. There is a justification too: if the features do not correlate with one another and the loss is quadratic, this formula is nothing other than the exact solution of the least squares problem. In practice this almost never happens, the assumption is completely impractical, but as a zeroth approximation — why not.
• Rough training on a small subsample. For example, for a linear classifier the direction vector can be roughly determined like this: take a random object of class $+1$ and a random object of class $-1$, and their difference as vectors in ${\mathbb{R}}^n$ gives the direction from the negative class to the positive one. If we take a few dozen such pairs and average them, we get a much better direction, roughly showing how one class differs from the other. A huge number of heuristics of this sort can be invented.
• Multistart — repeated runs from different random initial approximations, choosing the best solution. Especially appropriate when parallel computing is available: we run the iterative process in parallel from a random set of initial approximations, after some number of iterations choose the best, and can multistart once more from the neighbourhood of that best one.

The order in which objects are presented

Stochastic gradient means that objects are chosen one at a time in random order. Must it be random, and can some non-randomness not be brought in to speed up convergence? Here too there is a huge number of heuristics.

• Shuffling — for example, alternating objects of different classes. The assumption is: if we take an object similar to others that are classified reliably, the gradient step will be small, so it is worth taking objects that differ substantially from each other, and objects of different classes do differ.
• Selecting objects by margin. One can look at the margins introduced at the beginning of the lecture. If the margin is large, the object can be skipped altogether: all is well with it anyway, and gradient steps on it are useless. And if the margin is large in absolute value and negative, things are so bad with it that they will not get better: by taking it into account we will only spoil the classifier. Noise is noise; it would be best ignored altogether. In general: the larger the loss ${\epsilon }_i$, the higher the probability of choosing $x_i$.

Around these ideas very many hybrid heuristic methods of speeding up convergence have been invented.

Choosing the gradient step

Another question I have so far swept under the carpet is the gradient step $h$. The theory of optimizing convex functions says that the step should be decreased as the iteration number grows, but neither too fast nor too slowly:

in particular, one can set $h_t=\frac{1}{t}$. One can, but should not, because it works badly; there are theorems about this for convex functions, and if you optimize something more complex, God knows how the step should be decreased. We would like automatic methods for choosing the gradient step.

The method of steepest gradient descent. Take the formula of one iteration of the gradient method and treat the step size $h$ as a one-dimensional parameter over which we minimize the loss:

Sometimes this problem can even be solved analytically: differentiate with respect to $h$, set the derivative to zero and find a formula for $h^{\ast }$. But this cannot be done for all loss functions, and then one has to use numerical methods of one-dimensional minimization: the Fibonacci method, bisection, or Newton's method if we can compute the first and second derivatives. The main thing is no fanaticism: there is no need to optimize over $h$ super-precisely; it is enough to find $h$ roughly for one iteration.

Second-order methods. One can try to improve the first-order method by bringing second derivatives into it. In the first-order method the gradient is multiplied by the step size, while in the Newton–Raphson method it also has to be multiplied on the left by the inverse Hessian — the matrix of second partial derivatives:

The real nightmare here is that the Hessian must not only be computed but also inverted, and it is an $N\times N$ matrix, where $N$ is the dimension of the parameter vector. Too laborious; nobody does it. How can this formula nevertheless be roughly approximated so as to keep the features of a second-order method, in particular the high rate of convergence?

The diagonal Levenberg–Marquardt method. The Hessian is often diagonally dominant. Let us use a trick: assume that there is nothing in it except the diagonal, that is, the mixed partial derivatives with respect to different variables $w_j$, $w_{j^{\prime }}$ are zeros. Of course that is not so, but we simply will not compute them. Inverting a diagonal matrix is a pleasure: it is just ones divided by the second derivatives. But here too not everything is smooth: the second derivatives may be zero or close to zero, they will end up in the denominator, and we will take a random jump due to computational errors. So for stability one more parameter $\mu$ is added to the denominator:

We wanted to get rid of the gradient step parameter, and another one has appeared. But the good news is that the learning rate $h$ in the Levenberg–Marquardt method can be set to one: one divided by the second derivative of the loss is in fact an almost optimal gradient step. And $\mu$ is needed where the second derivatives are close to zero — at inflection points, where the loss in the space of $w$ is close to linear. Second-order methods do not work there: a second-order method is a parabolic approximation of the criterion being optimized, and it works well where we have almost reached the minimum; that is where the second derivatives give very fast convergence. This is how the gradient step is usually chosen, sometimes by simpler methods, for example by estimating the second derivatives with finite differences; many techniques have been invented.

Advantages and drawbacks of the SG method

So the stochastic gradient method seems wonderful, simple and good, but without additional effort — speeding up convergence, choosing the step, choosing the initial approximation — it may not work very well: converge slowly, get stuck in local extrema, and many other unpleasant things can happen to it. That is why over seventy years it has been surrounded by heuristics.

Its advantages:

1. Simplicity: it is relatively easy to implement.
2. Universality: it suits any models $a(x,w)$ and any loss functions $\mathcal{L}(w,x)$ — as long as the loss can be differentiated with respect to the parameter vector and is smooth. How exactly it is built does not much matter to the method: it is very flexible and can solve a very wide class of problems, not only linear models, both classification and regression.
3. Streaming: it can learn from a stream of data.
4. Suits big data: it can reach a decent solution having processed only a small part of the training sample, and it often turns out faster and better than more complex and resource-hungry second-order methods.

A few decades ago it was thought in optimization that first-order methods are not very interesting: they converge slowly compared with second-order methods. In machine learning it turned out otherwise. When big data appeared — the functional is a sum over the objects of the sample, and there are millions of objects, sometimes billions — stochastic gradient, especially with momentum, the accelerated gradient and other techniques, turned out to be a real treasure, and one practically does not have to bother with second-order methods. What remains of them is essentially the methods of choosing the gradient step, and real second-order methods can be used only in spaces of small dimension, where the Hessian can be computed and inverted.

The drawback, perhaps, is that one still has to put together a set of heuristics, and that is an art: do not forget overfitting, getting stuck, divergence. True, there are libraries for training neural networks with this very stochastic gradient method inside, and so much invented over the last decades is built into it that we no longer even look inside. We regard it as a reliable engine, a black box that can optimize functions of millions and billions of terms. And this account was so that you have a rough idea of what is under the hood of this black box.

Multicollinearity and regularization

The stochastic gradient method is surrounded by heuristics, but never mind that. Here, though, is a fundamental problem that shows up in linear models — the problem of multicollinearity; for linear models it is also the problem of overfitting. The point is this: some features may be linearly dependent. In particular, if there are fewer objects than features, there will certainly be linearly dependent features.

What is bad about that? Suppose we have built a binary classifier $a(x,w)=\mathrm{sign}⟨w,x⟩$. Multicollinearity is the condition that there is a vector $u\in {\mathbb{R}}^n$ such that for any object $x$ the scalar product $⟨u,x⟩=0$: with some coefficients $u$ the features can be added up so that the result is always zero. Then the solution is not unique, and therefore unstable. One can always add $\gamma u$ to $w$: in the scalar product with $x$ the addition collapses, $⟨u,x⟩=0$, and the classifier will work exactly as the original one.

The scalar $\gamma$ can be anything; it can be sent to plus or minus infinity, and the weights of the features will run off to plus or minus infinity. But the features are collinear only on the training data. Introduce small errors into their measurement, and the addition $\gamma ⟨u,x⟩$ instantly stops being zero. Then the scalar product flies off to an unpredictable place. Our learning method was looking for the optimum of the functional under multicollinearity and landed at some unknown point of the linear manifold of solutions. The optimum is unstable. This shows up as follows:

• the coefficients go haywire: weights $|w_j|$ that are too large in absolute value and of different signs appear;
• the discriminant function is unstable: change the vector $x$ ever so slightly, and its value flies off who knows where;
• overfitting: $Q(X^{\ell })\ll Q(X^k)$. On the training sample the classifier is fitted and everything seems fine, but on an independent control sample everything is bad because of this instability.

How to fight it? In the course we will see at least three ways:

• regularization: $\Vert w\Vert \to \min$ (weight decay);
• feature selection: $f_1,\dots ,f_n\to f_{j_1},\dots ,f_{j_m}$, $m\ll n$;
• feature transformation: $f_1,\dots ,f_n\to g_1,\dots ,g_m$, $m\ll n$.

Let us now discuss the simplest. As doctors say, the treatment is symptomatic. If there are infinitely many solutions, a whole manifold, and many of them run off to plus or minus infinity, let us choose among all the solutions the most stable one in the sense of the minimum norm of the vector $w$: minimize the coordinates of $w$ in absolute value as far as possible. That is why the technique is called weight decay. The functional leads us into a linear manifold of solutions, and we say: from this manifold we will choose the solution closest to the origin. A second criterion has appeared, the squared norm of the vector $w$, and we add the two criteria with some coefficient $\tau$.

This technique is also called stabilization, but most often regularization. It is completely standard, and for linear models simply obligatory, because we do not know whether there will be overfitting due to multicollinearity. Let us see what happens to the gradient vector: a small correction $\tau w$ has been added to it — because it is a squared norm, and with $\frac{1}{2}$ at that:

The modification of the gradient step in the SG method is, one could say, an edit to one line of code:

We simply multiply the vector $w$ from the current iteration by a number less than one — there it is, decay, shrinking, reducing the norm of the vector $w$ — and take the gradient step in the same direction as before. A very simple implementation.

The probabilistic statement of the classification problem

Another fundamentally important paradigm in machine learning, and an important message of today's lecture, is the probabilistic statement of the problem. Suppose the Cartesian product $X\times Y$ is a probability space with some density, and the «object — answer» pairs, both from the training sample and from the control one, are generated by this joint distribution independently: i.i.d. (independent identically distributed). In mathematical statistics such samples are called simple. Our task is to estimate the joint density from the sample: if we have an estimate of it, we will be able to compute the probability of the class $y$ for any object $x$.

What does it mean to find the density $p(x,y)$? We solve optimization problems over parameters all the time, so parameters must be introduced into the density. The density factorizes by the formula of conditional probability, which is always true, and we introduce the parameter $w$ into the model of the conditional probability of the class — that is convenient for us. The alternative approach, where the parametrization goes the other way round, we will also consider later. How to specify the model itself we will put off for now; first let us deal with the optimization criterion.

The maximum likelihood method is one of the fundamental methods of mathematical statistics and machine learning, and we will use it far and wide as soon as we talk about a probabilistic model of data generation. Once upon a time there was a density, objects were generated from it, and now we observe them. The direct problem is to generate objects from the density; we have to solve the inverse one: from the observed objects say what density could have generated them.

The likelihood of the data is the density of the whole sample with our data substituted into it. Why is it a product of densities over the individual objects? Because we assumed independent generation: the density of an independently generated sample equals the product of the densities of the individual pairs. We substitute the factorization of the density, and we are entitled to drop the factor $p(x_i)$: it does not depend on $w$. We need the density in order to say for which $w$ from the parametric family the generated sample is most likely. This is the maximum likelihood principle: choose the parameter of the generating density for which the observed data is most likely.

A product over all the objects of the sample is inconvenient for optimization, so one takes its logarithm: the maximum point of a functional and of its logarithm is the same. This is a completely general technique. As soon as you read about a probabilistic model of data generation, expect optimization of the log-likelihood two lines later. We get a sum over all the objects of the sample, and the logarithm of the model of the conditional class probability is nothing other than a loss function — more precisely, minus the logarithm: here we maximize, while losses are usually minimized.

Likelihood and the approximated empirical risk

There is a close link between maximizing likelihood and minimizing empirical risk — I have in fact already let it slip. Let us look at the same classification problem in two ways. On the one hand, we can build a probabilistic model of classification, for simplicity with two classes. On the other hand, we can say: we have no probabilistic hypothesis, no probabilistic nature of the data. We build a separating surface between the classes, we take a purely geometric view of the problem, and as the model we fix a parametric family of discriminant functions. Let us write out the two criteria that arise in these two paradigms.

Maximizing likelihood, where $P(y\mid x,w)$ is the model of the conditional probability of the class $y$:

Minimizing the approximated empirical risk, where $g(x,w)$ is the model of the separating surface, $Y=\{\pm 1\}$:

Both are sums over all the objects of the sample, where some loss is written for each object, and it is tempting to say that they are the same thing. Each term depends on only one object, so they can be equated independently.

Then the probabilistic approach is merely the construction of a separating surface with such a loss function. A dual view appears. People of different faiths may come: some believe in probability and that the world around us is probabilistic, while others regard probability as a strange thing — where in life have you seen algebras of measurable events? — and use pure geometry. Both are right, and they can be reconciled. You do not believe in probabilities, you have written out a separating surface $g$ and a loss function $L$? Look at this equality: I will take its exponential and say that this is the probabilistic model of the data. And from the other side a probabilist comes and says: «I am modelling the data». And we tell them: «Nothing of the kind, you are building a separating surface with such and such a loss function». These are two legitimate, equal, dual views of what we are doing. Whether to use the probabilistic paradigm or not is a matter of taste, a matter of which model of the world we believe in.

The probabilistic meaning of regularization

Since with this transition we can freely jump between the probabilistic and non-probabilistic paradigms, let us ask: what, in probabilistic language, is the regularization that arose for linear models?

Suppose there is a probabilistic data model $P(y\mid x,w)$. Who is $w$, where does it come from? Until now $w$ was a non-random parameter that we found from the sample by maximum likelihood. Now suppose that $w$ is also a random vector, generated from its own distribution. This is a two-level model of data generation. The space $W$ is, as it were, the space of the practical problems we face. The Lord God generates a problem from a prior distribution and slips it to us, and then a training sample is generated from this problem; we look at the data and try to determine the model. If we know the model of problem generation the Lord God uses, we can bring something extra into the solution.

The hyperparameters $\gamma$ are again a non-random vector, parameters over parameters. Since the likelihood has been formalized this way, nothing prevents us from writing out the maximum likelihood principle, which in this case is called the maximum a posteriori principle; it too is used very often in probabilistic statements of machine learning. The logarithm of the joint density gives two terms. The first is our old acquaintance, the probabilistic loss (if taken with a minus). The second, the logarithm of the prior distribution, does not depend on the data; it is a criterion on the parameters. That is, it is a regularizer.

Again two paradigms. In one we say: we have peeked into the heavenly chancellery and know from which distribution problems are slipped to us for practical solution. In the other, we are solving an ill-posed problem and have therefore added a reasonable restraint, a regularizer. These are two equivalent languages, and which of them to formalize the problem statement in is a matter of taste.

The Gaussian and Laplace priors

Since we have agreed that prior distributions, or rather their logarithms, can play the role of regularizers, let us make these assumptions: the weights $w_j$ are independent, the expectation of each is zero, and the variance is some constant $C$:

We take the logarithm and — an astonishing thing — get our acquaintance, the $L_2$ regularizer of the linear model. On the other hand, in probabilistic language we say that the parameter vector is somewhere in the neighbourhood of zero and does not run too far from it: there it is, a probabilistic assumption, a constraint on the variance. The regularization coefficient turns out to be inverse to the variance: the larger the variance, the weaker the regularization, and vice versa. A small variance means we press the weights hard towards zero, so the regularization coefficient is large.

$C$ is a hyperparameter, $\tau =\frac{1}{C}$ is the regularization coefficient.

The Laplace distribution resembles the Gaussian, but instead of the square it has the absolute norm — the sum of absolute values. This is a very interesting regularizer. $L_2$ regularization, weight decay, shrinks the coefficients: they decrease but do not reach zero. With $L_1$ regularization the effect of zeroing coefficients arises: the corresponding feature is simply not used; the regularizer switches features off. We will have many occasions to talk about this in the course.

We will discuss logistic regression and the example of the credit scoring problem another time.

Seminar