Machine Learning 1. Lecture 7. Generalizations of Linear Models

What this lecture is about

Today's lecture continues the previous one. Let us recall what we did last time: we considered the regression problem, the linear regression model and least squares, and within these restrictions learned a lot of interesting things about the solution. Now we begin to drop these restrictions and extend the problem statement to a wider class of problems. Today we will have nonlinear generalizations of regression; situations where linear models can be used but least squares has to be abandoned — that happens too; and non-standard, non-quadratic loss functions. In short, all sorts of generalizations — and the nasty things that happen when we deviate from the classics told in the last lecture.

The lecture will consist of a sequence of relatively separate topics, but they will have a lot in common. One idea will run through the whole lecture: a complex problem can be reduced to an iterative process, at each step of which a simpler problem is solved — as a rule, that very problem of multiple linear regression that we already know how to solve.

Nonlinear regression

Let us start straight with the problem statement. The objects are described by numerical features, that is, they are vectors in ${\mathbb{R}}^n$; the answers are real numbers, that is, we recover a numerical dependence of $y$ on $x$. An unknown regression dependence, we assume, exists, but we do not know it. Unlike in the last lecture, the regression model is an arbitrary parametric family of functions $f(x,\alpha )$, nonlinear.

Let us see how far one can get in solving the problem without fixing $f$: let it be some arbitrary nonlinear model. We still want to use least squares — minimize over $\alpha$ the sum of squared deviations of the model's answer from the correct value $y_i$. Last time the model was linear: we differentiated, set the derivative to zero, got a system of equations and studied its properties. That will not work here: in general the problem cannot be solved analytically, only by numerical methods.

We already know one numerical method — the first-order, gradient method, in particular the stochastic gradient method. It can simply be taken and applied to this problem. Full stop, lecture over, let us go home. But no such luck: let us try to apply a second-order method.

The Newton–Raphson method

About Newton methods, in particular the Newton–Raphson method, it is known that they converge much faster, and we can probably profit from that somehow and learn to solve the problem faster than first-order methods allow. Why the formula of the method is exactly so I will not tell: that is from the course on numerical optimization methods.

The main thing is that the iterative process consists of steps in the direction of the antigradient. The gradient points in the direction of fastest increase of the function, and we are solving a minimization problem, so we take it with a minus. But unlike first-order gradient methods, the gradient vector has to be multiplied by the inverse Hessian. The Hessian is the matrix of second partial derivatives: if the gradient is a vector of dimension $p$, equal to the dimension of the parameter vector $\alpha$, then the Hessian is a $p\times p$ matrix. The gradient step $h_t$ also appears, but in the Newton–Raphson method it is known that the step can be of the order of one; it can even be set to one. At least some relief: there is no need to choose the step size. On the other hand, technical difficulties arise with the Hessian, and even more with inverting it. Let us gradually work out how to simplify this problem.

First let us write out the components of the gradient vector — the partial derivatives of the functional with respect to ${\alpha }_j$. Note: unlike in the stochastic gradient method, I do not get rid of the summation sign; we must honestly carry out the sum over all the objects of the sample; there will be no small steps on each object.

We still keep the function $f$ in general form and want to carry it through to the end, to arrive at an iterative process into which any $f$ can be substituted to get an optimization algorithm. Let us differentiate once more, with respect to ${\alpha }_k$, and get the components of the Hessian:

The expression consists of two terms. With the first we can do nothing, and the second, it is claimed, can be neglected. Why? Let us use a technique physicists love when they meet nonlinear models: linearization. If a function is smooth enough, in the neighbourhood of some point it can be approximately replaced by a linear one. We are at the point ${\alpha }^t$, where $t$ is the iteration number, and compute the Hessian at this point. Let us, in a small neighbourhood of ${\alpha }^t$, replace the nonlinear model by its linear approximation — a Taylor expansion:

The value of the function at the point ${\alpha }^t$ is a known constant; it can be computed. Then the linear term: ${\alpha }_j$ is a variable, ${\alpha }_j^t$ its fixed value at the point ${\alpha }^t$, and the derivative is also taken at a fixed point. The terms of second, third and higher orders we discard. If we substitute this linearization into the second term of the Hessian, the first derivative with respect to a parameter will turn out to be a constant, and the second derivative zero. That is, under linearization the second term vanishes. And the first does not change: the first derivatives at the point ${\alpha }^t$ of the original function and of its linearization are the same.

The Gauss–Newton method

Let us look carefully at the first term that remains. What does it remind us of? A product of two matrices — the matrices of partial derivatives of the model $f$ with respect to all the parameters ${\alpha }_j$. Let us move to matrix form.

The matrix of first derivatives has size $\ell \times p$: the rows correspond to objects, the columns to parameters. It is like the «objects–features» matrix, only the role of the features is played by the partial derivatives of the model with respect to each of the parameters. $f_t$ is the vector of values of the model on all the objects of the training sample. Let us rewrite the formula of one iteration of the Newton–Raphson method with linearization: the second big term with second derivatives we have discarded, and what remains is in effect the product $F_t^{\mathsf{T}}{F}_t$ — it is what we invert; then we multiply by the gradient. We get an expression called the Gauss–Newton method. It is no longer the original Newton–Raphson method: we have replaced the Hessian by an approximation of it. But this method too is known to be of second order, and the rate of convergence has not suffered from linearization.

And now look what we have got. It is the pseudoinverse matrix from the last lecture! There it was built from the «objects–features» matrix as the result of solving the least squares problem for a linear regression model. Let us do some reverse engineering: what linear regression problem have we solved? Clearly, we approximated the vector of differences «model minus correct answers» by a linear combination of the vectors of derivatives of the model with respect to all its parameters. In meaning this also resembles a gradient step: we try to approximate the difference between the model and the vector of correct answers by an optimally chosen sum of the model's derivatives with respect to the parameters. So we have reduced the solution of the nonlinear regression problem to an iterative process, at each step of which a multiple linear regression problem of its own is solved. And that we already know how to solve, and all the methods of the last lecture — the singular value decomposition, principal component analysis if the dimension needs reducing, regularization — fit here and can be applied.

Logistic regression

What we have just considered is nonlinear regression in which least squares remained, and that led us to a second-order method. In the second topic we will keep the model linear but take a linear classification model. The loss function will of course be different: there is no point in using least squares when solving a classification problem. Let us look at logistic regression. We have already talked about first-order methods, about stochastic gradient, and they solve logistic regression well; today we will consider an alternative path — again the Newton–Raphson iterative process.

A reminder: a two-class classification problem; we denote the classes $-1$ and $+1$. The linear classification model is the sign of a scalar product; $w$ is the direction vector of the separating hyperplane, which points towards class $+1$, $x$ are vectors in ${\mathbb{R}}^n$, and $w$ is a vector of the same dimension. The number of errors on the training sample is written through the notion of the margin: the margin shows the signed distance from the $i$th point of the training sample to the separating surface; if there is no error, the sign is positive, if there is, negative. We replace the threshold loss with an upper bound. Many different such bounds can be made: SVM is one bound, and if we take the smooth function $\ln (1+e^{-M})$, we get a method called logistic regression. The Newton–Raphson method here is an alternative at least because, first, the loss function is different, and second, we are not putting regularization in here yet, though nothing prevents doing so.

The Newton–Raphson method for logistic regression

Again a second-order optimization method, and again I do not remove the sum, no stochastic gradient: we optimize the sum over all the objects of the sample.

Let us differentiate honestly. The elements of the gradient, the vector of first derivatives $\nabla Q(w^t)$:

The elements of the Hessian, the matrix of second derivatives $Q^{\prime \prime }(w^t)$:

where ${\sigma }_i=\sigma (y_i{w}^{\mathsf{T}}{x}_i)$ and $\sigma (z)=\frac{1}{1+e^{-z}}$ is the sigmoid function.

What is ${\sigma }_i$? It is the value of the sigmoid at the margin of the $i$th object. This quantity has a probabilistic meaning, and we will talk about it now, and somewhere in the middle of the lecture I will prove why it is so. Let us look at the Hessian: up to the factor $(1-{\sigma }_i){\sigma }_i$ it is the componentwise notation of the product of two matrices $F^{\mathsf{T}}F$. That is, we again recognize the matrix we have to invert: the Hessian again takes on a familiar form.

Again a reduction to multiple linear regression

Let us move to matrix notation. $F$ is the «objects–features» matrix: the rows are the $\ell$ objects of the training sample, the columns are the features. From the expressions $(1-{\sigma }_i){\sigma }_i$ let us form a diagonal matrix $D$: these values are on the diagonal, all other elements are zero. The factors $(1-{\sigma }_i){\sigma }_i$ get in the way in the product of two matrices, and inserting them there amounts to putting the diagonal matrix $D$ in the middle, between $F^{\mathsf{T}}$ and $F$.

Our aggregate, which gives an iteration of the Newton–Raphson method, is again a construction very similar to the pseudoinverse matrix, only with $D$ in the middle. To come to the form $(F^{\mathsf{T}}F{)}^{-1}{F}^{\mathsf{T}}$, we have to get rid of $D$ somehow, and that is not hard: take the root of $D$ — a diagonal matrix with the corresponding roots on its diagonal — and multiply $F$ by it on the left and $F^{\mathsf{T}}$ on the right. Transforming a little, we will again see that at each iteration we solve a multiple linear regression problem, with the elements of $D$ turning out to be the weights of the objects in the least squares statement. The regression model here is linear, and the answers $y_i=\pm 1$ are modified: they are divided by ${\sigma }_i$.

An amazing construction. Looking at it, one could try the following: take $y_i=\pm 1$, approximate them by a linear model and solve the classification problem by ordinary least squares for linear regression. In principle one can, nobody forbids it, though it will not be a very good solution: fitting a linear regression to the numbers $-1$ and $+1$, to class labels. But here we are told that the class labels must be normalized — divided by the probabilities of correct classification.

If an object is classified wrongly, ${\sigma }_i$ is small, and the modified answer $y_i/{\sigma }_i$ runs off from zero to plus or minus infinity. In this way we tune more strongly to the classification of such an object: if the model errs on the object, we pay more attention to it by increasing its $|y_i|$. For classification it does not matter what the absolute value of the model's answer turns out to be: we take its sign, and only the sign interests us. So we can play with the size of the $y_i$ that the regression is tuned to, in order to solve a classification problem.

But that is not all: the modification has also affected the weights of the objects. If ${\sigma }_i$ is the probability of correct classification, which objects will have the largest weight? Those where ${\sigma }_i=\frac{1}{2}$. These are the objects on the boundary between the classes, with the most uncertain classification. So, building such a linear regression model, we also pay more attention to the objects lying on the boundary between the classes. We are more attentive to the objects that have turned out to be hard for the current linear classification model at the current iteration. And hardness is of two types: either we err on the object, or the object is borderline and we are unsure of its classification. It is mainly on such objects that we train the model at the next iteration.

Iteration after iteration, a linear regression problem is solved at each step. Hence, by the way, the name: why is logistic regression a regression if it solves a classification problem? That is why. The second-order method I have told you about is the textbook, classical method for solving the logistic regression problem. And it is also a regression because it predicts the probabilities of objects belonging to a class, and a probability is no longer a discrete quantity but a continuous one. This too brings logistic regression closer to regression, although a classification problem is being solved.

The IRLS method

Let us put everything together and write the pseudocode of the iterative process, so as not to get lost in the intermediate computations.

The zeroth approximation is that very unfortunate way of using regression for classification: we took our $y_i=\pm 1$ and substituted them into the ordinary solution of the least squares problem for multiple linear regression, with the pseudoinverse matrix from the last lecture. It is a rough solution of the classification problem by regression, but it will do as a zeroth approximation. Then at each iteration, for the current approximation of the vector $w$, we compute ${\sigma }_i$ — the probabilities of correctly classifying the objects of the training sample — compute the diagonal matrix, multiply the features by it on the left and reweight the answers, multiplying them by a factor that depends only on ${\sigma }_i$. So we have reweighted both the objects and the answers. The gradient step in the Newton–Raphson method can be set to one, or it can be finely tuned at each iteration if one wishes. What happens at the eighth step is the solution of a multiple linear regression problem. We have reduced a new, harder problem to the successive solution of problems we already know how to solve, as was done in the last lecture.

The generalized additive model

The third topic falls somewhat out of the main narrative, but it is an important thing, and it had to be told somewhere; now is the best moment. We want to generalize the linear regression problem, but in a limited version: not to introduce an arbitrary nonlinear dependence on $\alpha$, but to say that all we lack for a linear model is to transform each individual feature a little — to bring it into a form in which a linear model can be built over the transformed features.

The best example is the applied problem I gave in the first lecture: predicting real estate prices. The price of a flat in Moscow is made up of many factors: floor area, the number of rooms, a balcony, and God knows what else. Some features are numerical, others binary, like having a balcony. But what about a feature such as the distance to the metro? Intuition suggests that the greater the distance to the metro, the smaller the contribution to the price: the closer to the metro, the higher the price of the flat. One could put this distance into a linear regression, and probably a negative coefficient would come out. But are we sure that the distance to the metro makes a literally additive, linear contribution? Rather not: it should be some decreasing function, and we would like to find out what it looks like — to train it too, tune it, optimize it. There are quite a lot of such problems: one wants first to perform a separate learnable nonlinear transformation on each feature separately, and then build a linear regression model over the transformed features.

The model is the sum over all the features of their transformations. A transformation may also be linear — simply multiplying the feature $f_j$ by the coefficient ${\alpha }_j$, and then it is a special case, the ordinary linear regression model. But we want to allow nonlinear transformations of individual features, and we need to train the functions ${\phi }_j$ — functions of a scalar argument with some parameter. Let us write least squares for such a model and train the ${\phi }_j$ in turn: training one, we fix all the others; having trained it, we take another feature and train its ${\phi }_j$, and so on until the process converges to something. The fixed other functions form a construction we will call $z_i$, and we get the problem of training a one-dimensional regression model, no longer on $y_i$ but on the differences between $y_i$ and the sum of all the other, fixed transformations of the features.

Regularization can also be added here. We may assume that the function ${\phi }_j$ transforming the feature $f_j$ is smooth enough, that is, does not differ too much from a linear one. That very decreasing function that converts the distance to the metro into a unit price — how much the price of a flat decreases with every hundred metres from the metro — is monotonically decreasing, probably not too complicated and not too different from a linear dependence. We do not mean it to have any humps, ridges or sinusoids. How do we introduce such an additional constraint — that we want a function as simple as possible in the sense of smoothness? An example is the regularizer that is the integral of the squared second derivative: for a linear function the second derivative is zero, and the integral is zero; a small non-smoothness, a small nonlinearity lead to small values of this penalty term.

There are methods that directly solve the one-dimensional problem with such a regularizer, splines for example. One can also do it differently — by kernel smoothing with the Nadaraya–Watson formula. Choosing the smoothing window width there is roughly the same as choosing the regularization coefficient $\tau$: the larger $\tau$, the smoother the function we get, and the wider the window, the smoother the function we get. So the one-dimensional problem is fairly easy to solve. Moreover, the one-dimensional problem is easy to visualize: one can plot the approximation of the one-dimensional dependence of $z_i$ on $f_j(x_i)$ and see whether the chosen function ${\phi }_j$ really describes it well.

The backfitting method

We get a method into which visual quality control of the solution can be inserted at each step. It is called backfitting, and it is rather old: Hastie and Tibshirani invented it back in 1986. It is good precisely in problems like predicting real estate prices, where one wants to check with one's eyes how meaningful the nonlinear transformations of each feature come out. For the distance to the metro you expect to get some nontrivial ${\phi }_j$ — a monotonically decreasing function, and it is interesting to look at it. But for the floor area of the flat you hardly expect any clever nonlinear transformation and may well leave the linear dependence as it is.

Look how interesting. In the previous two topics we also built an iterative process, and at each iteration the nonlinear problem was replaced by the classical linear regression problem. Here too a replacement takes place, but a different one: solving the problem at each iteration reduces to solving a one-dimensional nonlinear regression. There we replaced a nonlinear model by a linear one at each iteration; here we replace a multidimensional problem by a one-dimensional one.

Let us go through the code. The initial approximation: again one can solve the multiple linear regression problem and say that ${\phi }_j$ is simply multiplying $f_j$ by the coefficient ${\alpha }_j$. This is the first approximation; we will substitute it when computing $z_i$ — what the one-dimensional model is tuned to for each feature; the target for us is now $z_i$, but only at this iteration; at each iteration it is different. I write the one-dimensional problem in general form, though it is implied that you choose yourselves: splines, Nadaraya–Watson kernel smoothing, some other approximation — anything. But in any of these methods a parameter that has to be optimized will appear one way or another. With too strong smoothing you will get a linear model, and there will be no nonlinearity in the feature transformation. With too weak regularization you may overfit: zero errors on the objects of the training sample, but a too complex, rapidly oscillating function $\phi$. So the regularization has to be chosen by looking at how the function behaves and how well it passes through the sample points. In this way, step by step, by successively complicating the solution — gradually decreasing $\tau$ — we build the model.

The probabilistic meaning of least squares

Now the most important part of today's lecture, in my view — generalized linear models. But first let us sort out the names. The backfitting method is a generalized additive model; do not confuse it with the generalized linear model, which we are about to discuss. The names are similar; the meaning is completely different. Let us start with the close link between least squares and maximum likelihood and understand what probabilistic meaning least squares carries. It is nontrivial.

We will build a model of the data, assuming that the $y_i$ are predicted by our model up to the noise ${\epsilon }_i$. This noise is the residual between the model and the true unknown target dependence. We will assume that the model is so good that the residuals have certain probabilistic properties: they can be considered normally distributed with zero expectation. Zero expectation is the assumption that the model is unbiased: it approximates the target dependence so that the errors are now positive, now negative, with no systematic bias. And another assumption is that we know the variances ${\sigma }_i^2$ or can estimate them. And in general the variances on different objects may differ: on some object we can predict more accurately, and on some the noise variance is larger, and we cannot claim high prediction accuracy there. Sometimes a stricter assumption is made, that all the variances ${\sigma }_i$ are equal — this is called the homoscedasticity assumption for the residuals. We are now in a more general statement: unknown and different variances are assigned to the residuals ${\epsilon }_i$.

An equivalent notation: the $y_i$ have a normal distribution with some expectations ${\mu }_i$ and unknown variances ${\sigma }_i^2$, and the expectation of $y_i$ is our model. Such is the probabilistic view of regression: regression describes the expectations of the $y_i$, up to noise, and about the distribution of this noise we assume Gaussianity.

Let us apply maximum likelihood, since we have ended up in a probabilistic statement. The random variables here are the ${\epsilon }_i$; let us write out their joint density: they are distributed with zero expectation and variances ${\sigma }_i^2$.

Where are the parameters $w$ here? ${\epsilon }_i=y_i-a(x_i,w)$: on the one hand the ${\epsilon }_i$ are random variables, and on the other they are functions of $w$, and the problem can be solved as an optimization problem over $w$. The most interesting thing: if we maximize the log-likelihood, we get a sum of logarithms, the logarithm and the exponential destroy each other, and we unexpectedly see that the maximum likelihood principle is equivalent to least squares. ${\epsilon }_i^2$ is the square of the difference between the model's answer and $y_i$. This is least squares, with a small generalization: we also ended up with weights of the objects.

And this too is an interesting insight: it turns out that in least squares the weight of an object is nothing other than one divided by the noise variance on this object. The larger the noise, the smaller the weight: if we understand that a good prediction is impossible on some objects, there is a large noise variance there — the weight of such an object in the least squares criterion should be reduced proportionally. So least squares is not taken out of thin air: behind it stands the implicit assumption that our model is so good that it is unbiased, has zero expected errors, and that the noise is Gaussian.

A question arises: if for some reason we are not sure that the errors are Gaussian, do we then have no right to use least squares? Generally speaking, yes: something else has to be used, having made another, more plausible assumption about the form of the error distribution. And one more consideration: if we predict an integer quantity rather than a real-valued one, can we expect the errors to be Gaussian? Most likely not, and we should understand what error distribution would be more adequate in this case. What should be used instead of least squares if the $y_i$ are not Gaussian, in particular if they are discrete-valued?

The generalized linear model

Let us generalize least squares so as to approximate dependences by models for which we are not sure that the errors are Gaussian and can write out some alternative formula for the error distribution. Such a generalization is called the generalized linear model: we do want to stay within linear models.

Let us start again from the normal linear model for the expectation: we assume that the $y_i$ are normally distributed and their expectation is described by a linear model. This is the standard statement, closely tied to least squares. Now we generalize: instead of the Gaussian distribution we take a very wide class of distributions called the exponential family. But everything has to be paid for. For widening the class of distributions we will pay with the fact that the parameter ${\theta }_i$, which we will model linearly, is now not equal to the expectation of $y_i$ but is tied to it through a known function — the link function.

The exponential family of distributions

The exponential family is specified up to two parameters, ${\theta }_i$ and ${\varphi }_i$, called the location and the dispersion respectively. They are some analogues of the expectation and the variance, but they are exactly these only for the Gaussian distribution; since not exactly, they have other names. And the exponential family is also specified by two functional parameters: the function $c(\theta )$ and the function $h(y,\varphi )$.

Look what a formula, from who knows where. The main thing is that it is the exponential of something, and inside it there is a form linear in $y_i$: $y_i$ is multiplied by ${\theta }_i$; the function $c$ depends on $\theta$ but not on $\varphi$ or $y$; and the function $h$ does not depend on $\theta$. A kind of decomposition has happened: the exponential of a sum is a product of exponentials, and the distribution can be represented as a product of two functions, one of which does not depend on $\theta$, and the other depends on $y$ only linearly. I will now give examples, and it will turn out that very many well-known distributions are special cases of the exponential family.

What is important for us to know about the exponential family? First, for a random variable generated by the exponential family one can compute the expectation and the variance. I omit the calculation: it is not hard; the corresponding integral has to be computed. It turns out that the expectation is $c^{\prime }$, the derivative of the function $c$ of the scalar parameter $\theta$. Hence $\theta$ is expressed as a function of the expectation: the inverse function of the derivative. We conclude that $c^{\prime }$ must be monotonic so that its inverse can be taken. This $g$ is the link function, linking the expectation with the parameter ${\theta }_i$. And ${\theta }_i$ we declared from the start to be modelled by our linear model — not the expectation, as before in the simple case of least squares. For the generalization we have to pay with the appearance of the function $g$. But we know how to compute it: as soon as a particular distribution from the exponential family is written out, we know the function $c$, can differentiate it and find the inverse. That will be the link function. So it does not create any special technical difficulties for us.

Examples of distributions from the exponential family

Let us consider four examples.

What has to be done to get everything under the exponential in the normal distribution? Take the logarithm of the normalizing factor and bring it inside, and expand the squared difference and regroup the brackets a little. $c^{\prime }(\theta )=\theta$ is the identity function, and such a miracle happens only in the Gaussian distribution as a special case of the exponential family. And ${\varphi }_i$ here really does turn out to equal the variance — hence the name «dispersion parameter»: in the Gaussian case it equals the variance, and for other distributions it may be something else.

Why is the Bernoulli distribution interesting? Because it is the two-class classification problem. When the residuals are normal, it is a regression problem, where we are entitled to apply least squares. But if we are solving a classification problem, we are not entitled to apply least squares: we cannot assume that the $y_i$ are Gaussian. They are distributed by Bernoulli: one class has probability ${\mu }_i$, the other $1-{\mu }_i$, that is, $y_i=1$ with probability ${\mu }_i$ and $y_i=0$ with probability $1-{\mu }_i$; the formula writes the same thing as a product. Bringing this notation to the exponential form is very easy: we artificially introduce an exponential — first put the expression under a logarithm, then take its exponential — and indeed $y_i$ turns out to be multiplied by something depending only on the expectation. That is ${\theta }_i$, and $c$ can be expressed through it. There is no dispersion parameter here; it equals one.

The binomial distribution is the number of successes in a Bernoulli scheme. It arises in problems where we want to approximate the probability of a discrete event, knowing that the maximum number of times the event could have happened is known and equals $n_i$. On each object there is some counter: the maximum number of times is known to us, and we observed the event some number of times, and these data are the $y_i$ to be predicted: how many times the event will happen on the given object. Since we predict an integer quantity — the number of successes in a Bernoulli scheme — it has a known distribution, the binomial. At first sight the formula even has the combinatorial coefficients $C_{n_i}^{y_i}$, but it is easy to show that this too is an exponential family: we put everything under the exponential and are surprised to find a linear term with exactly the same formula as for the Bernoulli distribution. Well, not quite the same: there it was one, and here it is a known integer $n_i$. The dispersion parameter does not need estimating either.

The Poisson distribution is a limiting special case of the previous one: we do not know how many times the event could have happened, but we predict some integer quantity. For example, the number of clicks on an advertising banner: CTR (click-through rate) in many problems of internet advertising, in problems related to recommender systems and in many other places, an integer quantity arises, and it is unknown what its upper limit is. Unlike the binomial distribution, where the maximum number of successes is known in advance, here it is unknown how many times users could have clicked, yet one wants to predict the probability of a click. Then the Poisson distribution may be more appropriate. Again we put everything under the exponential, find $y_i$ multiplied by something depending only on the expectation — and that gives the link function of the expectation ${\mu }_i$ with ${\theta }_i$, which we approximate by a linear model. And ${\varphi }_i$ again equals one. That is, in three discrete distributions very often needed in practice, the dispersion parameter does not need estimating, and that is a pleasant simplification.

I have shown four examples of how one proves that a known distribution is a special case of the exponential family: put everything under the exponential, collect like terms, and it becomes obvious what the functions $c$ and $h$ equal. But this cannot be done with every distribution.

The list of useful distributions that are not exponential is in fact much longer, and it is rather an exception for a distribution to be exponential. But in practice it is exponential distributions that arise very often. Once again: this is about the distribution of the quantity $y$ that we are going to predict. If some theory tells us that the predicted quantity is distributed according to some distribution — ${\chi }^2$, Dirichlet, beta, gamma, Poisson and so on — this theoretical knowledge should be treated carefully. Since we are predicting a quantity that is not Gaussian, we have no right to use least squares: we must use maximum likelihood explicitly, modelling the distribution of the $y_i$.

Maximizing the likelihood for GLM

Let us do this in the general case, for an arbitrary distribution from the exponential family, and see along the way why it has exactly this form.

where ${\theta }_i$ depends linearly on $w$: ${\theta }_i=⟨x_i,w⟩={\sum }_{j=1}^n{w}_j{f}_j(x_i)$.

We substitute the formula of the exponential distribution and take the logarithm at once: the logarithm and the exponential destroy each other. But I immediately used one more trick. The function $h$ depends on $y$ and $\varphi$; if ${\varphi }_i$ is a constant, one for example, then $h$ can be thrown out; it will not depend on the parameters. Moreover, if we assume homoscedasticity — the dispersions are the same for all objects, that is, there is no reason in the model to think that on some objects the errors will be much larger than on others — the formula simplifies: only the first terms remain, and the terms with $h$ go away, because they depend neither on the model parameters nor on the objects of the sample. This cannot always be done; look carefully: if this term does not drop out, it has to be taken into account carefully. And once again let me remind you where $w$ is here: ${\theta }_i$ is the linear model itself, the sum of all the features of the object with the weights $w_j$, so the maximization is over $w$, on which $\theta$ depends.

A good, pleasant functional, and it has to be optimized. We have already done this trick twice today; let us do it a third time — the Newton–Raphson method:

You can already guess what will happen next. We honestly differentiate with respect to $w_j$ and $w_k$.

The components of the gradient vector $\nabla Q(w)$:

The components of the Hessian $Q^{\prime \prime }(w)$:

And again we will have to invert a matrix very similar to the product $F^{\mathsf{T}}F$, only a factor has crept in that depends on the object but not on the features $j$ and $k$. Again this will be a diagonal matrix inserted between $F^{\mathsf{T}}$ and $F$, and we will probably now see iteratively reweighted least squares again.

From the «objects–features» matrix we form a new, reweighted one, and we reweight the answers: subtract the expectation $c^{\prime }({\theta }_i)$ and divide by something like the spread — the root of the variance $\mathsf{D}{y}_i={\varphi }_i{c}^{\prime \prime }({\theta }_i)$; clearly this is the modified vector of answers. And the Newton–Raphson method again leads to iterative reweighting, that is, to solving a linear regression problem. I will not interpret it further: that is somewhat more complicated than for logistic regression, but the meaning is roughly the same. Again we have reduced the problem to a series of multiple linear regression problems. Let us write the pseudocode once more, now for the more general statement.

The essence of the process is the same. The zeroth approximation is ordinary least squares, as if it could be applied; it will do. Then we compute the linear model, the weights, reweight; the gradient step $h_t$ can be set to one or chosen. At the eighth step is the pseudoinverse matrix, that is, the solution of the multiple linear regression problem for the answers and the «objects–features» matrix transformed by these formulas. That is all — very similar to what we got for logistic regression. And this is no accident: substituting $c^{\prime }$ and $c^{\prime \prime }$ of the Bernoulli distribution, we get exactly IRLS for logistic regression. Logistic regression is a special case of the generalized linear model.

Logistic regression as a special case of GLM

Let us consider logistic regression as a special case of GLM and understand the main thing: where the sigmoid function in logistic regression comes from, and where the probabilistic interpretation comes from. It turns out that the only source of this is the probabilistic assumption that we predict $y_i$ distributed by Bernoulli. The essence of this assumption is just one thing: we assumed that for each object there exists a probability of class $+1$. We made no other probabilistic assumptions. And this probability is also the parameter of the Bernoulli distribution and also the expectation. Everything else follows automatically from the general theory.

${\theta }_i$ is our linear model, ${\mu }_i$ is the expectation, that is, the probability of class $+1$ for the $i$th object. The link between them is the link function; if we write the inverse function — wow, it is the sigmoid function! So that is where it came from. And if we write the maximum likelihood functional, that very logarithmic loss function pops out — a direct consequence of the link function, that is, of a single probabilistic assumption: $y_i$ has a probability of belonging to a class. From it both the sigmoid and the logarithmic loss in the functional equivalent to maximum likelihood popped out at once. That is, these are not heuristics but consequences of the probabilistic assumption that a probability of class $+1$ exists for each object of the training sample. Look what a harmless probabilistic assumption we made, and what interesting consequences we got in the form of a particular loss function. This, in my view, is a little bit fantastic, a very beautiful mathematical consequence. And at last you understand where logistic regression came from.

Let us do the calculation I described in words. Here we have to do a little alchemy with the coding of the classes: in the Bernoulli distribution $y_i\in \{0,1\}$, while in linear models we are used to denoting the classes $-1$ and $+1$. One can be recoded into the other, for example, by the formula ${\tilde{y}}_i=2y_i-1$. For the object $x_i$ the probability of belonging to class $+1$ is the expectation, through the link function it is $\sigma ({\theta }_i)$. And the probability of class $0$ is $1-{\mu }_i=\sigma (-{\theta }_i)$, thanks to the property of the sigmoid function $1-\sigma (z)=\sigma (-z)$. From these two equalities we conclude that the probability of ${\tilde{y}}_i$ for $x_i$ is the sigmoid of the margin: $P({\tilde{y}}_i\mid x_i)=\sigma (⟨w,x_i⟩{\tilde{y}}_i)$. So that is what it is: the margin, written through ${\tilde{y}}_i$ taking the values $\pm 1$.

The maximum of likelihood — which is also the minimum of the log-loss criterion — can be written in two equivalent ways. On the one hand, it is the sum over all the objects of the sample of the logarithms of the probability of the class $y_i$ for the object $x_i$, written through the Bernoulli distribution. This form of the criterion is often called log-loss, and in very many places, in many binary classification problems, everything starts by writing out this criterion. The other way of writing is through the logarithmic loss of the margin. One can write it either way when one meets a probabilistic statement of the problem. This, in my view, is the central insight of today's lecture about logistic regression; now one can say that you know the main things about it. And the optimization here can be done by the second-order method, IRLS, or by the stochastic gradient method I have already shown you.

One more small addition: the value the linear model gives on an object also has an interesting probabilistic meaning. It is ${\theta }_i$, which is also $\ln \frac{{\mu }_i}{1-{\mu }_i}$, that is, the logarithm of the odds ratio — the logarithm of the ratio of the probability of the object being in class $+1$ to the probability of it being in the other class. It turns out that the numerical value of the scalar product is not only the distance from the object to the separating hyperplane but also the logarithm of the odds ratio. A very interesting probabilistic interpretation of the linear model. And all this, let me remind you, follows from just a single assumption: the classes $y_i$ have a Bernoulli distribution, and the Bernoulli distribution has a single parameter, which is also the expectation, the probability of class $+1$.

When least squares will not do

The last topic is other cases where we have no right to use least squares, although we would very much like to: it is such a developed method, all the ready numerical methods exist for it.

Let us start with the fact that we do not always want to minimize squared errors. Least squares has a solution for the case where the model is a constant; we have already discussed it, including when deriving the Nadaraya–Watson formula. The model is just a scalar parameter $\alpha$, with no dependence on $x$; differentiating with respect to $\alpha$ and setting the derivative to zero, we immediately get that $\alpha$ is the mean value of $y_i$ on the objects of the training sample. And if we want to minimize not the sum of squared deviations but the sum of absolute values? Solving the analogous problem for a constant, we will see that the solution is the median of the $y_i$, not the arithmetic mean.

What is good about the median? It is much more robust to outliers. If a very large error has crept in among the $y_i$ — for example, instead of $y_i$ something a million times larger was taken — the arithmetic mean will jerk very strongly in that direction: it is sensitive to large errors in measuring the $y_i$. But the median does not care; it will not jerk anywhere: the anomalously large value will end up last in the order statistics, and in the middle, most likely, the same value or its neighbour, very close to it, will remain. Hence the property: the method of least absolute values, with such a loss function, is preferable in problems where extremely large outliers are expected. But what does that mean? That the distribution of the ${\epsilon }_i$ is not Gaussian. Assuming the presence of outliers, we assume the Laplace distribution — with heavy tails, also exponential, but depending on the absolute value of the error rather than its square.

Quantile regression

A generalization: what if the loss function is not just the absolute value but a skewed absolute value, and the penalties for positive errors are smaller than the penalties for negative ones?

For example, in the problem of forecasting sales volumes. If we made an underforecast — the goods were not delivered — we created a shortage and did not get the trade markup: those are large losses. And if an overforecast — too much was delivered, funds were frozen, the goods sat on the shelf but still sell — the loss is smaller. Such situations arise in different applied problems. If we want to measure losses literally in roubles, they will be a linear, not a quadratic, function of the forecast errors: a square of roubles has no economic meaning, accountants tell us. Why use least squares in such problems? Skewed absolute values should be used.

What will the constant solution be? Instead of the median a quantile of the distribution will appear — hence the name quantile regression. If we take the order statistics of the values $y_i$ and in them take not the middle value but the $\gamma$th, where $\gamma$ is computed from the slopes of the two branches of the loss function, then this solution is optimal. For a primitive, constant model the solution of quantile regression is a quantile, that is, an element of the order statistics of the values of $y$. The more heavily an underforecast is penalized, the higher up the order statistics the optimal forecast.

And if the regression model is linear, quantile regression is easily written out so that it reduces to a linear programming problem — minimizing a linear criterion under linear equality and inequality constraints. For example, like this: represent the residual $y_i-⟨w,x_i⟩$ as the difference of two non-negative variables, the positive and negative parts, as we did for LASSO,

The problem of forecasting sales volumes is in practice quite often solved precisely with quantile regression and linear programming. It can be even trickier; I believe I showed this in the first lecture when telling about forecasting sales volumes: the loss function may be piecewise linear made not of two pieces but, for example, of five. With different sizes of overforecast the loss factors change: the more we overbought, the more new factors come into play, and the losses begin to grow with a different coefficient of proportionality. And the large losses from an underforecast are tied to the loss of the trade markup and of customer loyalty.

Robust loss functions

There are also loss functions invented in order to ignore large errors in measuring the $y_i$. For example, instead of the quadratic loss we want to take a function with horizontal asymptotes. Many such functions can be invented.

Here, as an example, is Meshalkin's function. The problem with such a loss function is also not least squares; it too is, one could say, maximization of likelihood, and one can even write out the error distribution — it will of course be not Gaussian but some other. But one can also not bother with probabilistic interpretations and simply say: we took such a loss function, and we solve the problem by the same Newton–Raphson method. I will not write it out for you a fourth time.

There is also Barron's parametric family of loss functions: by changing the parameter one can get both the parabolic loss (at the value 2) and loss functions with horizontal asymptotes. Within one parametrization one can change the model's attitude to noise in measuring the $y_i$. In general robust regression is a big research topic; thick books on robust regression have been written, and I have a separate lecture on robust models; for now only a spoiler, and the lecture will be in another semester. Many functions have been invented — Meshalkin's, Huber's, Tukey's, Andrews' and many other surnames — and most of them have the property that their branches, unlike the parabola, are either horizontal or have a linear order of growth. That is, they are more tolerant of large errors and, when the model is tuned, simply pay no attention to them.

SVM regression I have also already told you about in the lecture on SVM: it too is an example where in a regression problem the classical loss function was replaced by a piecewise linear one and all the properties of the SVM method were obtained — support objects and a quadratic programming problem; everything is very similar to SVM for classification. So we have run through the cases where something other than least squares has to be used, and as a rule this follows from the understanding that our errors are not Gaussian.

Summary

We have come to the point where an iterative process has appeared three times, no, even four times: in nonlinear regression, in logistic regression, in the generalized linear model and in problems with non-quadratic loss functions. That is today's main message: complex problems can be reduced to an iterative process — to a series of solutions of simpler linear problems that we already know how to solve.

• Logistic regression is not regression but classification; the Newton–Raphson method leads to IRLS.
• The generalized linear model (GLM) powerfully generalizes ordinary and logistic regression; the Newton–Raphson method leads to IRLS.
• The generalized additive model (GAM, backfitting) reduces to a series of one-dimensional smoothings.
• The probabilistic meaning of the loss function: it is $-\ln$ of the probabilistic model of the data.

Today's main topic is still the generalized linear model. It is applied in very many places, both in machine learning theory and in practice, and is well studied; I gave you a simplified account. Thick surveys and books have been written about it, but for a general understanding of what was today, this is enough for you.

Seminar