master

CS229 Spring 2022 2 5 Kernel methods 48 5.1 Feature maps . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.2 LMS (least mean squares) with features . . . . . . . . . . . . . 49 5.3 LMS with the kernel trick . . . . . . . . . . . . . . . . . . . . 49 5.4 Properties of kernels . . . . . . . . . . . . . . . . . . . . . . . 53 6 Support vector machines 59 6.1 Margins: intuition . . . . . . . . . . . . . . . . . . . . . . . . . 59 6.2 Notation (option reading) . . . . . . . . . . . . . . . . . . . . 61 6.3 Functional and geometric margins (option reading) . . . . . . 61 6.4 The optimal margin classifier (option reading) . . . . . . . . . 63 6.5 Lagrange duality (optional reading) . . . . . . . . . . . . . . . 65 6.6 Optimal margin classifiers: the dual form (option reading) . . 68 6.7 Regularization and the non-separable case (optional reading) . 72 6.8 The SMO algorithm (optional reading) . . . . . . . . . . . . . 73 6.8.1 Coordinate ascent . . . . . . . . . . . . . . . . . . . . . 74 6.8.2 SMO . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 II Deep learning 79 7 Deep learning 80 7.1 Supervised learning with non-linear models . . . . . . . . . . . 80 7.2 Neural networks . . . . . . . . . . . . . . . . . . . . . . . . . . 82 7.3 Backpropagation . . . . . . . . . . . . . . . . . . . . . . . . . 91 7.3.1 Preliminary: chain rule . . . . . . . . . . . . . . . . . . 92 7.3.2 One-neuron neural networks . . . . . . . . . . . . . . . 92 7.3.3 Two-layer neural networks: a low-level unpacked com- putation . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.3.4 Two-layer neural network with vector notation . . . . . 96 7.3.5 Multi-layer neural networks . . . . . . . . . . . . . . . 98 7.4 Vectorization over training examples . . . . . . . . . . . . . . 98 III Generalization and regularization 101 8 Generalization 102 8.1 Bias-variance tradeoff . . . . . . . . . . . . . . . . . . . . . . . 104 8.1.1 A mathematical decomposition (for regression) . . . . . 109 8.2 The double descent phenomenon . . . . . . . . . . . . . . . . . 110

5 Let’s start by talking about a few examples of supervised learning prob- lems. Suppose we have a dataset giving the living areas and prices of 47 houses from Portland, Oregon: Living area (feet 2 ) Price (1000 $ s) 2104 400 1600 330 2400 369 1416 232 3000 540 . . . . . . We can plot this data: 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 100 200 300 400 500 600 700 800 900 1000 housing prices square feet price (in $1000) Given data like this, how can we learn to predict the prices of other houses in Portland, as a function of the size of their living areas? To establish notation for future use, we’ll use x ( i ) to denote the “input” variables (living area in this example), also called input features , and y ( i ) to denote the “output” or target variable that we are trying to predict (price). A pair ( x ( i ) , y ( i ) ) is called a training example , and the dataset that we’ll be using to learn—a list of n training examples { ( x ( i ) , y ( i ) ); i = 1 , . . . , n } —is called a training set . Note that the superscript “( i )” in the notation is simply an index into the training set, and has nothing to do with exponentiation. We will also use X denote the space of input values, and Y the space of output values. In this example, X = Y = R . To describe the supervised learning problem slightly more formally, our goal is, given a training set, to learn a function h : X 7→ Y so that h ( x ) is a “good” predictor for the corresponding value of y . For historical reasons, this

8 confusion, we will drop the θ subscript in h θ ( x ), and write it more simply as h ( x ). To simplify our notation, we also introduce the convention of letting x 0 = 1 (this is the intercept term ), so that h ( x ) = d X i =0 θ i x i = θ T x, where on the right-hand side above we are viewing θ and x both as vectors, and here d is the number of input variables (not counting x 0 ). Now, given a training set, how do we pick, or learn, the parameters θ ? One reasonable method seems to be to make h ( x ) close to y , at least for the training examples we have. To formalize this, we will define a function that measures, for each value of the θ ’s, how close the h ( x ( i ) )’s are to the corresponding y ( i ) ’s. We define the cost function : J ( θ ) = 1 2 n X i =1 ( h θ ( x ( i ) ) - y ( i ) ) 2 . If you’ve seen linear regression before, you may recognize this as the familiar least-squares cost function that gives rise to the ordinary least squares regression model. Whether or not you have seen it previously, let’s keep going, and we’ll eventually show this to be a special case of a much broader family of algorithms. 1.1 LMS algorithm We want to choose θ so as to minimize J ( θ ). To do so, let’s use a search algorithm that starts with some “initial guess” for θ , and that repeatedly changes θ to make J ( θ ) smaller, until hopefully we converge to a value of θ that minimizes J ( θ ). Specifically, let’s consider the gradient descent algorithm, which starts with some initial θ , and repeatedly performs the update: θ j := θ j - α ∂ ∂θ j J ( θ ) . (This update is simultaneously performed for all values of j = 0 , . . . , d .) Here, α is called the learning rate . This is a very natural algorithm that repeatedly takes a step in the direction of steepest decrease of J . In order to implement this algorithm, we have to work out what is the partial derivative term on the right hand side. Let’s first work it out for the

11 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 0 100 200 300 400 500 600 700 800 900 1000 housing prices square feet price (in $1000) If the number of bedrooms were included as one of the input features as well, we get θ 0 = 89 . 60 , θ 1 = 0 . 1392, θ 2 = - 8 . 738. The above results were obtained with batch gradient descent. There is an alternative to batch gradient descent that also works very well. Consider the following algorithm: Loop { for i = 1 to n , { θ j := θ j + α ( y ( i ) - h θ ( x ( i ) ) ) x ( i ) j , (for every j ) (1.2) } } By grouping the updates of the coordinates into an update of the vector θ , we can rewrite update (1.2) in a slightly more succinct way: θ := θ + α ( y ( i ) - h θ ( x ( i ) ) ) x ( i ) In this algorithm, we repeatedly run through the training set, and each time we encounter a training example, we update the parameters according to the gradient of the error with respect to that single training example only. This algorithm is called stochastic gradient descent (also incremental gradient descent ). Whereas batch gradient descent has to scan through the entire training set before taking a single step—a costly operation if n is large—stochastic gradient descent can start making progress right away, and

14 Finally, to minimize J , let’s find its derivatives with respect to θ . Hence, ∇ θ J ( θ ) = ∇ θ 1 2 ( Xθ - ~ y ) T ( Xθ - ~ y ) = 1 2 ∇ θ ( ( Xθ ) T Xθ - ( Xθ ) T ~ y - ~ y T ( Xθ ) + ~ y T ~ y ) = 1 2 ∇ θ ( θ T ( X T X ) θ - ~ y T ( Xθ ) - ~ y T ( Xθ ) ) = 1 2 ∇ θ ( θ T ( X T X ) θ - 2( X T ~ y ) T θ ) = 1 2 ( 2 X T Xθ - 2 X T ~ y ) = X T Xθ - X T ~ y In the third step, we used the fact that a T b = b T a , and in the fifth step used the facts ∇ x b T x = b and ∇ x x T Ax = 2 Ax for symmetric matrix A (for more details, see Section 4.3 of “Linear Algebra Review and Reference”). To minimize J , we set its derivatives to zero, and obtain the normal equations : X T Xθ = X T ~ y Thus, the value of θ that minimizes J ( θ ) is given in closed form by the equation θ = ( X T X ) - 1 X T ~ y. 3 1.3 Probabilistic interpretation When faced with a regression problem, why might linear regression, and specifically why might the least-squares cost function J , be a reasonable choice? In this section, we will give a set of probabilistic assumptions, under which least-squares regression is derived as a very natural algorithm. Let us assume that the target variables and the inputs are related via the equation y ( i ) = θ T x ( i ) + ( i ) , 3 Note that in the above step, we are implicitly assuming that X T X is an invertible matrix. This can be checked before calculating the inverse. If either the number of linearly independent examples is fewer than the number of features, or if the features are not linearly independent, then X T X will not be invertible. Even in such cases, it is possible to “fix” the situation with additional techniques, which we skip here for the sake of simplicty.

17 0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x y 0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x y 0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 x y Instead, if we had added an extra feature x 2 , and fit y = θ 0 + θ 1 x + θ 2 x 2 , then we obtain a slightly better fit to the data. (See middle figure) Naively, it might seem that the more features we add, the better. However, there is also a danger in adding too many features: The rightmost figure is the result of fitting a 5-th order polynomial y = ∑ 5 j =0 θ j x j . We see that even though the fitted curve passes through the data perfectly, we would not expect this to be a very good predictor of, say, housing prices ( y ) for different living areas ( x ). Without formally defining what these terms mean, we’ll say the figure on the left shows an instance of underfitting —in which the data clearly shows structure not captured by the model—and the figure on the right is an example of overfitting . (Later in this class, when we talk about learning theory we’ll formalize some of these notions, and also define more carefully just what it means for a hypothesis to be good or bad.) As discussed previously, and as shown in the example above, the choice of features is important to ensuring good performance of a learning algorithm. (When we talk about model selection, we’ll also see algorithms for automat- ically choosing a good set of features.) In this section, let us briefly talk about the locally weighted linear regression (LWR) algorithm which, assum- ing there is sufficient training data, makes the choice of features less critical. This treatment will be brief, since you’ll get a chance to explore some of the properties of the LWR algorithm yourself in the homework. In the original linear regression algorithm, to make a prediction at a query point x (i.e., to evaluate h ( x )), we would: 1. Fit θ to minimize ∑ i ( y ( i ) - θ T x ( i ) ) 2 . 2. Output θ T x . In contrast, the locally weighted linear regression algorithm does the fol- lowing: 1. Fit θ to minimize ∑ i w ( i ) ( y ( i ) - θ T x ( i ) ) 2 . 2. Output θ T x .

20 is called the logistic function or the sigmoid function . Here is a plot showing g ( z ): -5 -4 -3 -2 -1 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z g(z) Notice that g ( z ) tends towards 1 as z → ∞ , and g ( z ) tends towards 0 as z → -∞ . Moreover, g(z), and hence also h ( x ), is always bounded between 0 and 1. As before, we are keeping the convention of letting x 0 = 1, so that θ T x = θ 0 + ∑ d j =1 θ j x j . For now, let’s take the choice of g as given. Other functions that smoothly increase from 0 to 1 can also be used, but for a couple of reasons that we’ll see later (when we talk about GLMs, and when we talk about generative learning algorithms), the choice of the logistic function is a fairly natural one. Before moving on, here’s a useful property of the derivative of the sigmoid function, which we write as g 0 : g 0 ( z ) = d dz 1 1 + e - z = 1 (1 + e - z ) 2 ( e - z ) = 1 (1 + e - z ) · 1 - 1 (1 + e - z ) = g ( z )(1 - g ( z )) . So, given the logistic regression model, how do we fit θ for it? Following how we saw least squares regression could be derived as the maximum like- lihood estimator under a set of assumptions, let’s endow our classification model with a set of probabilistic assumptions, and then fit the parameters via maximum likelihood.

21 Let us assume that P ( y = 1 | x ; θ ) = h θ ( x ) P ( y = 0 | x ; θ ) = 1 - h θ ( x ) Note that this can be written more compactly as p ( y | x ; θ ) = ( h θ ( x )) y (1 - h θ ( x )) 1 - y Assuming that the n training examples were generated independently, we can then write down the likelihood of the parameters as L ( θ ) = p ( ~ y | X ; θ ) = n Y i =1 p ( y ( i ) | x ( i ) ; θ ) = n Y i =1 ( h θ ( x ( i ) ) ) y ( i ) ( 1 - h θ ( x ( i ) ) ) 1 - y ( i ) As before, it will be easier to maximize the log likelihood: ‘ ( θ ) = log L ( θ ) = n X i =1 y ( i ) log h ( x ( i ) ) + (1 - y ( i ) ) log(1 - h ( x ( i ) )) How do we maximize the likelihood? Similar to our derivation in the case of linear regression, we can use gradient ascent. Written in vectorial notation, our updates will therefore be given by θ := θ + α ∇ θ ‘ ( θ ). (Note the positive rather than negative sign in the update formula, since we’re maximizing, rather than minimizing, a function now.) Let’s start by working with just one training example ( x, y ), and take derivatives to derive the stochastic gradient ascent rule: ∂ ∂θ j ‘ ( θ ) = y 1 g ( θ T x ) - (1 - y ) 1 1 - g ( θ T x ) ∂ ∂θ j g ( θ T x ) = y 1 g ( θ T x ) - (1 - y ) 1 1 - g ( θ T x ) g ( θ T x )(1 - g ( θ T x )) ∂ ∂θ j θ T x = ( y (1 - g ( θ T x )) - (1 - y ) g ( θ T x ) ) x j = ( y - h θ ( x )) x j

22 Above, we used the fact that g 0 ( z ) = g ( z )(1 - g ( z )). This therefore gives us the stochastic gradient ascent rule θ j := θ j + α ( y ( i ) - h θ ( x ( i ) ) ) x ( i ) j If we compare this to the LMS update rule, we see that it looks identical; but this is not the same algorithm, because h θ ( x ( i ) ) is now defined as a non-linear function of θ T x ( i ) . Nonetheless, it’s a little surprising that we end up with the same update rule for a rather different algorithm and learning problem. Is this coincidence, or is there a deeper reason behind this? We’ll answer this when we get to GLM models. 2.2 Digression: the perceptron learning algo- rithn We now digress to talk briefly about an algorithm that’s of some historical interest, and that we will also return to later when we talk about learning theory. Consider modifying the logistic regression method to “force” it to output values that are either 0 or 1 or exactly. To do so, it seems natural to change the definition of g to be the threshold function: g ( z ) = 1 if z ≥ 0 0 if z < 0 If we then let h θ ( x ) = g ( θ T x ) as before but using this modified definition of g , and if we use the update rule θ j := θ j + α ( y ( i ) - h θ ( x ( i ) ) ) x ( i ) j . then we have the perceptron learning algorithn . In the 1960s, this “perceptron” was argued to be a rough model for how individual neurons in the brain work. Given how simple the algorithm is, it will also provide a starting point for our analysis when we talk about learning theory later in this class. Note however that even though the perceptron may be cosmetically similar to the other algorithms we talked about, it is actually a very different type of algorithm than logistic regression and least squares linear regression; in particular, it is difficult to endow the perceptron’s predic- tions with meaningful probabilistic interpretations, or derive the perceptron as a maximum likelihood estimation algorithm.

23 1 1.5 2 2.5 3 3.5 4 4.5 5 -10 0 10 20 30 40 50 60 x f(x) 1 1.5 2 2.5 3 3.5 4 4.5 5 -10 0 10 20 30 40 50 60 x f(x) 1 1.5 2 2.5 3 3.5 4 4.5 5 -10 0 10 20 30 40 50 60 x f(x) 2.3 Another algorithm for maximizing ‘ ( θ ) Returning to logistic regression with g ( z ) being the sigmoid function, let’s now talk about a different algorithm for maximizing ‘ ( θ ). To get us started, let’s consider Newton’s method for finding a zero of a function. Specifically, suppose we have some function f : R 7→ R , and we wish to find a value of θ so that f ( θ ) = 0. Here, θ ∈ R is a real number. Newton’s method performs the following update: θ := θ - f ( θ ) f 0 ( θ ) . This method has a natural interpretation in which we can think of it as approximating the function f via a linear function that is tangent to f at the current guess θ , solving for where that linear function equals to zero, and letting the next guess for θ be where that linear function is zero. Here’s a picture of the Newton’s method in action: In the leftmost figure, we see the function f plotted along with the line y = 0. We’re trying to find θ so that f ( θ ) = 0; the value of θ that achieves this is about 1.3. Suppose we initialized the algorithm with θ = 4 . 5. Newton’s method then fits a straight line tangent to f at θ = 4 . 5, and solves for the where that line evaluates to 0. (Middle figure.) This give us the next guess for θ , which is about 2.8. The rightmost figure shows the result of running one more iteration, which the updates θ to about 1.8. After a few more iterations, we rapidly approach θ = 1 . 3. Newton’s method gives a way of getting to f ( θ ) = 0. What if we want to use it to maximize some function ‘ ? The maxima of ‘ correspond to points where its first derivative ‘ 0 ( θ ) is zero. So, by letting f ( θ ) = ‘ 0 ( θ ), we can use the same algorithm to maximize ‘ , and we obtain update rule: θ := θ - ‘ 0 ( θ ) ‘ 00 ( θ ) . (Something to think about: How would this change if we wanted to use Newton’s method to minimize rather than maximize a function?)

24 Lastly, in our logistic regression setting, θ is vector-valued, so we need to generalize Newton’s method to this setting. The generalization of Newton’s method to this multidimensional setting (also called the Newton-Raphson method) is given by θ := θ - H - 1 ∇ θ ‘ ( θ ) . Here, ∇ θ ‘ ( θ ) is, as usual, the vector of partial derivatives of ‘ ( θ ) with respect to the θ i ’s; and H is an d -by- d matrix (actually, d +1 - by - d+1, assuming that we include the intercept term) called the Hessian , whose entries are given by H ij = ∂ 2 ‘ ( θ ) ∂θ i ∂θ j . Newton’s method typically enjoys faster convergence than (batch) gra- dient descent, and requires many fewer iterations to get very close to the minimum. One iteration of Newton’s can, however, be more expensive than one iteration of gradient descent, since it requires finding and inverting an d -by- d Hessian; but so long as d is not too large, it is usually much faster overall. When Newton’s method is applied to maximize the logistic regres- sion log likelihood function ‘ ( θ ), the resulting method is also called Fisher scoring .

Chapter 3 Generalized linear models So far, we’ve seen a regression example, and a classification example. In the regression example, we had y | x ; θ ∼ N ( μ, σ 2 ), and in the classification one, y | x ; θ ∼ Bernoulli( φ ), for some appropriate definitions of μ and φ as functions of x and θ . In this section, we will show that both of these methods are special cases of a broader family of models, called Generalized Linear Models (GLMs). 1 We will also show how other models in the GLM family can be derived and applied to other classification and regression problems. 3.1 The exponential family To work our way up to GLMs, we will begin by defining exponential family distributions. We say that a class of distributions is in the exponential family if it can be written in the form p ( y ; η ) = b ( y ) exp( η T T ( y ) - a ( η )) (3.1) Here, η is called the natural parameter (also called the canonical param- eter ) of the distribution; T ( y ) is the sufficient statistic (for the distribu- tions we consider, it will often be the case that T ( y ) = y ); and a ( η ) is the log partition function . The quantity e - a ( η ) essentially plays the role of a nor- malization constant, that makes sure the distribution p ( y ; η ) sums/integrates over y to 1. A fixed choice of T , a and b defines a family (or set) of distributions that is parameterized by η ; as we vary η , we then get different distributions within this family. 1 The presentation of the material in this section takes inspiration from Michael I. Jordan, Learning in graphical models (unpublished book draft), and also McCullagh and Nelder, Generalized Linear Models (2nd ed.) . 25

26 We now show that the Bernoulli and the Gaussian distributions are ex- amples of exponential family distributions. The Bernoulli distribution with mean φ , written Bernoulli( φ ), specifies a distribution over y ∈ { 0 , 1 } , so that p ( y = 1; φ ) = φ ; p ( y = 0; φ ) = 1 - φ . As we vary φ , we obtain Bernoulli distributions with different means. We now show that this class of Bernoulli distributions, ones obtained by varying φ , is in the exponential family; i.e., that there is a choice of T , a and b so that Equation (3.1) becomes exactly the class of Bernoulli distributions. We write the Bernoulli distribution as: p ( y ; φ ) = φ y (1 - φ ) 1 - y = exp( y log φ + (1 - y ) log(1 - φ )) = exp log φ 1 - φ y + log(1 - φ ) . Thus, the natural parameter is given by η = log( φ/ (1 - φ )). Interestingly, if we invert this definition for η by solving for φ in terms of η , we obtain φ = 1 / (1 + e - η ). This is the familiar sigmoid function! This will come up again when we derive logistic regression as a GLM. To complete the formulation of the Bernoulli distribution as an exponential family distribution, we also have T ( y ) = y a ( η ) = - log(1 - φ ) = log(1 + e η ) b ( y ) = 1 This shows that the Bernoulli distribution can be written in the form of Equation (3.1), using an appropriate choice of T , a and b . Let’s now move on to consider the Gaussian distribution. Recall that, when deriving linear regression, the value of σ 2 had no effect on our final choice of θ and h θ ( x ). Thus, we can choose an arbitrary value for σ 2 without changing anything. To simplify the derivation below, let’s set σ 2 = 1. 2 We 2 If we leave σ 2 as a variable, the Gaussian distribution can also be shown to be in the exponential family, where η ∈ R 2 is now a 2-dimension vector that depends on both μ and σ . For the purposes of GLMs, however, the σ 2 parameter can also be treated by considering a more general definition of the exponential family: p ( y ; η, τ ) = b ( a, τ ) exp(( η T T ( y ) - a ( η )) /c ( τ )). Here, τ is called the dispersion parameter , and for the Gaussian, c ( τ ) = σ 2 ; but given our simplification above, we won’t need the more general definition for the examples we will consider here.

27 then have: p ( y ; μ ) = 1 √ 2 π exp - 1 2 ( y - μ ) 2 = 1 √ 2 π exp - 1 2 y 2 · exp μy - 1 2 μ 2 Thus, we see that the Gaussian is in the exponential family, with η = μ T ( y ) = y a ( η ) = μ 2 / 2 = η 2 / 2 b ( y ) = (1 / √ 2 π ) exp( - y 2 / 2) . There’re many other distributions that are members of the exponen- tial family: The multinomial (which we’ll see later), the Poisson (for mod- elling count-data; also see the problem set); the gamma and the exponen- tial (for modelling continuous, non-negative random variables, such as time- intervals); the beta and the Dirichlet (for distributions over probabilities); and many more. In the next section, we will describe a general “recipe” for constructing models in which y (given x and θ ) comes from any of these distributions. 3.2 Constructing GLMs Suppose you would like to build a model to estimate the number y of cus- tomers arriving in your store (or number of page-views on your website) in any given hour, based on certain features x such as store promotions, recent advertising, weather, day-of-week, etc. We know that the Poisson distribu- tion usually gives a good model for numbers of visitors. Knowing this, how can we come up with a model for our problem? Fortunately, the Poisson is an exponential family distribution, so we can apply a Generalized Linear Model (GLM). In this section, we will we will describe a method for constructing GLM models for problems such as these. More generally, consider a classification or regression problem where we would like to predict the value of some random variable y as a function of x . To derive a GLM for this problem, we will make the following three assumptions about the conditional distribution of y given x and about our model:

28 1. y | x ; θ ∼ ExponentialFamily( η ). I.e., given x and θ , the distribution of y follows some exponential family distribution, with parameter η . 2. Given x , our goal is to predict the expected value of T ( y ) given x . In most of our examples, we will have T ( y ) = y , so this means we would like the prediction h ( x ) output by our learned hypothesis h to satisfy h ( x ) = E[ y | x ]. (Note that this assumption is satisfied in the choices for h θ ( x ) for both logistic regression and linear regression. For instance, in logistic regression, we had h θ ( x ) = p ( y = 1 | x ; θ ) = 0 · p ( y = 0 | x ; θ ) + 1 · p ( y = 1 | x ; θ ) = E[ y | x ; θ ].) 3. The natural parameter η and the inputs x are related linearly: η = θ T x . (Or, if η is vector-valued, then η i = θ T i x .) The third of these assumptions might seem the least well justified of the above, and it might be better thought of as a “design choice” in our recipe for designing GLMs, rather than as an assumption per se. These three assumptions/design choices will allow us to derive a very elegant class of learning algorithms, namely GLMs, that have many desirable properties such as ease of learning. Furthermore, the resulting models are often very effective for modelling different types of distributions over y ; for example, we will shortly show that both logistic regression and ordinary least squares can both be derived as GLMs. 3.2.1 Ordinary least squares To show that ordinary least squares is a special case of the GLM family of models, consider the setting where the target variable y (also called the response variable in GLM terminology) is continuous, and we model the conditional distribution of y given x as a Gaussian N ( μ, σ 2 ). (Here, μ may depend x .) So, we let the ExponentialFamily ( η ) distribution above be the Gaussian distribution. As we saw previously, in the formulation of the Gaussian as an exponential family distribution, we had μ = η . So, we have h θ ( x ) = E [ y | x ; θ ] = μ = η = θ T x. The first equality follows from Assumption 2, above; the second equality follows from the fact that y | x ; θ ∼ N ( μ, σ 2 ), and so its expected value is given

29 by μ ; the third equality follows from Assumption 1 (and our earlier derivation showing that μ = η in the formulation of the Gaussian as an exponential family distribution); and the last equality follows from Assumption 3. 3.2.2 Logistic regression We now consider logistic regression. Here we are interested in binary classifi- cation, so y ∈ { 0 , 1 } . Given that y is binary-valued, it therefore seems natural to choose the Bernoulli family of distributions to model the conditional dis- tribution of y given x . In our formulation of the Bernoulli distribution as an exponential family distribution, we had φ = 1 / (1 + e - η ). Furthermore, note that if y | x ; θ ∼ Bernoulli( φ ), then E[ y | x ; θ ] = φ . So, following a similar derivation as the one for ordinary least squares, we get: h θ ( x ) = E [ y | x ; θ ] = φ = 1 / (1 + e - η ) = 1 / (1 + e - θ T x ) So, this gives us hypothesis functions of the form h θ ( x ) = 1 / (1 + e - θ T x ). If you are previously wondering how we came up with the form of the logistic function 1 / (1 + e - z ), this gives one answer: Once we assume that y condi- tioned on x is Bernoulli, it arises as a consequence of the definition of GLMs and exponential family distributions. To introduce a little more terminology, the function g giving the distri- bution’s mean as a function of the natural parameter ( g ( η ) = E[ T ( y ); η ]) is called the canonical response function . Its inverse, g - 1 , is called the canonical link function . Thus, the canonical response function for the Gaussian family is just the identify function; and the canonical response function for the Bernoulli is the logistic function. 3 3.2.3 Softmax regression Let’s look at one more example of a GLM. Consider a classification problem in which the response variable y can take on any one of k values, so y ∈ { 1 , 2 , . . . , k } . For example, rather than classifying email into the two classes 3 Many texts use g to denote the link function, and g - 1 to denote the response function; but the notation we’re using here, inherited from the early machine learning literature, will be more consistent with the notation used in the rest of the class.

30 spam or not-spam—which would have been a binary classification problem— we might want to classify it into three classes, such as spam, personal mail, and work-related mail. The response variable is still discrete, but can now take on more than two values. We will thus model it as distributed according to a multinomial distribution. Let’s derive a GLM for modelling this type of multinomial data. To do so, we will begin by expressing the multinomial as an exponential family distribution. To parameterize a multinomial over k possible outcomes, one could use k parameters φ 1 , . . . , φ k specifying the probability of each of the outcomes. However, these parameters would be redundant, or more formally, they would not be independent (since knowing any k - 1 of the φ i ’s uniquely determines the last one, as they must satisfy ∑ k i =1 φ i = 1). So, we will instead pa- rameterize the multinomial with only k - 1 parameters, φ 1 , . . . , φ k - 1 , where φ i = p ( y = i ; φ ), and p ( y = k ; φ ) = 1 - ∑ k - 1 i =1 φ i . For notational convenience, we will also let φ k = 1 - ∑ k - 1 i =1 φ i , but we should keep in mind that this is not a parameter, and that it is fully specified by φ 1 , . . . , φ k - 1 . To express the multinomial as an exponential family distribution, we will define T ( y ) ∈ R k - 1 as follows: T (1) =        1 0 0 . . . 0        , T (2) =        0 1 0 . . . 0        , T (3) =        0 0 1 . . . 0        , · · · , T ( k - 1) =        0 0 0 . . . 1        , T ( k ) =        0 0 0 . . . 0        , Unlike our previous examples, here we do not have T ( y ) = y ; also, T ( y ) is now a k - 1 dimensional vector, rather than a real number. We will write ( T ( y )) i to denote the i -th element of the vector T ( y ). We introduce one more very useful piece of notation. An indicator func- tion 1 {·} takes on a value of 1 if its argument is true, and 0 otherwise (1 { True } = 1, 1 { False } = 0). For example, 1 { 2 = 3 } = 0, and 1 { 3 = 5 - 2 } = 1. So, we can also write the relationship between T ( y ) and y as ( T ( y )) i = 1 { y = i } . (Before you continue reading, please make sure you un- derstand why this is true!) Further, we have that E[( T ( y )) i ] = P ( y = i ) = φ i . We are now ready to show that the multinomial is a member of the

31 exponential family. We have: p ( y ; φ ) = φ 1 { y =1 } 1 φ 1 { y =2 } 2 · · · φ 1 { y = k } k = φ 1 { y =1 } 1 φ 1 { y =2 } 2 · · · φ 1 - ∑ k - 1 i =1 1 { y = i } k = φ ( T ( y )) 1 1 φ ( T ( y )) 2 2 · · · φ 1 - ∑ k - 1 i =1 ( T ( y )) i k = exp(( T ( y )) 1 log( φ 1 ) + ( T ( y )) 2 log( φ 2 ) + · · · + 1 - ∑ k - 1 i =1 ( T ( y )) i log( φ k )) = exp(( T ( y )) 1 log( φ 1 /φ k ) + ( T ( y )) 2 log( φ 2 /φ k ) + · · · + ( T ( y )) k - 1 log( φ k - 1 /φ k ) + log( φ k )) = b ( y ) exp( η T T ( y ) - a ( η )) where η =      log( φ 1 /φ k ) log( φ 2 /φ k ) . . . log( φ k - 1 /φ k )      , a ( η ) = - log( φ k ) b ( y ) = 1 . This completes our formulation of the multinomial as an exponential family distribution. The link function is given (for i = 1 , . . . , k ) by η i = log φ i φ k . For convenience, we have also defined η k = log( φ k /φ k ) = 0. To invert the link function and derive the response function, we therefore have that e η i = φ i φ k φ k e η i = φ i (3.2) φ k k X i =1 e η i = k X i =1 φ i = 1 This implies that φ k = 1 / ∑ k i =1 e η i , which can be substituted back into Equa- tion (3.2) to give the response function φ i = e η i ∑ k j =1 e η j

32 This function mapping from the η ’s to the φ ’s is called the softmax function. To complete our model, we use Assumption 3, given earlier, that the η i ’s are linearly related to the x ’s. So, have η i = θ T i x (for i = 1 , . . . , k - 1), where θ 1 , . . . , θ k - 1 ∈ R d +1 are the parameters of our model. For notational convenience, we can also define θ k = 0, so that η k = θ T k x = 0, as given previously. Hence, our model assumes that the conditional distribution of y given x is given by p ( y = i | x ; θ ) = φ i = e η i ∑ k j =1 e η j = e θ T i x ∑ k j =1 e θ T j x (3.3) This model, which applies to classification problems where y ∈ { 1 , . . . , k } , is called softmax regression . It is a generalization of logistic regression. Our hypothesis will output h θ ( x ) = E[ T ( y ) | x ; θ ] = E      1 { y = 1 } 1 { y = 2 } . . . 1 { y = k - 1 } x ; θ      =      φ 1 φ 2 . . . φ k - 1      =         exp( θ T 1 x ) ∑ k j =1 exp( θ T j x ) exp( θ T 2 x ) ∑ k j =1 exp( θ T j x ) . . . exp( θ T k - 1 x ) ∑ k j =1 exp( θ T j x )         . In other words, our hypothesis will output the estimated probability that p ( y = i | x ; θ ), for every value of i = 1 , . . . , k . (Even though h θ ( x ) as defined above is only k - 1 dimensional, clearly p ( y = k | x ; θ ) can be obtained as 1 - ∑ k - 1 i =1 φ i .)

33 Lastly, let’s discuss parameter fitting. Similar to our original derivation of ordinary least squares and logistic regression, if we have a training set of n examples { ( x ( i ) , y ( i ) ); i = 1 , . . . , n } and would like to learn the parameters θ i of this model, we would begin by writing down the log-likelihood ‘ ( θ ) = n X i =1 log p ( y ( i ) | x ( i ) ; θ ) = n X i =1 log k Y l =1 e θ T l x ( i ) ∑ k j =1 e θ T j x ( i ) ! 1 { y ( i ) = l } To obtain the second line above, we used the definition for p ( y | x ; θ ) given in Equation (3.3). We can now obtain the maximum likelihood estimate of the parameters by maximizing ‘ ( θ ) in terms of θ , using a method such as gradient ascent or Newton’s method.

Chapter 4 Generative learning algorithms So far, we’ve mainly been talking about learning algorithms that model p ( y | x ; θ ), the conditional distribution of y given x . For instance, logistic regression modeled p ( y | x ; θ ) as h θ ( x ) = g ( θ T x ) where g is the sigmoid func- tion. In these notes, we’ll talk about a different type of learning algorithm. Consider a classification problem in which we want to learn to distinguish between elephants ( y = 1) and dogs ( y = 0), based on some features of an animal. Given a training set, an algorithm like logistic regression or the perceptron algorithm (basically) tries to find a straight line—that is, a decision boundary—that separates the elephants and dogs. Then, to classify a new animal as either an elephant or a dog, it checks on which side of the decision boundary it falls, and makes its prediction accordingly. Here’s a different approach. First, looking at elephants, we can build a model of what elephants look like. Then, looking at dogs, we can build a separate model of what dogs look like. Finally, to classify a new animal, we can match the new animal against the elephant model, and match it against the dog model, to see whether the new animal looks more like the elephants or more like the dogs we had seen in the training set. Algorithms that try to learn p ( y | x ) directly (such as logistic regression), or algorithms that try to learn mappings directly from the space of inputs X to the labels { 0 , 1 } , (such as the perceptron algorithm) are called discrim- inative learning algorithms. Here, we’ll talk about algorithms that instead try to model p ( x | y ) (and p ( y )). These algorithms are called generative learning algorithms. For instance, if y indicates whether an example is a dog (0) or an elephant (1), then p ( x | y = 0) models the distribution of dogs’ features, and p ( x | y = 1) models the distribution of elephants’ features. After modeling p ( y ) (called the class priors ) and p ( x | y ), our algorithm 34

35 can then use Bayes rule to derive the posterior distribution on y given x : p ( y | x ) = p ( x | y ) p ( y ) p ( x ) . Here, the denominator is given by p ( x ) = p ( x | y = 1) p ( y = 1) + p ( x | y = 0) p ( y = 0) (you should be able to verify that this is true from the standard properties of probabilities), and thus can also be expressed in terms of the quantities p ( x | y ) and p ( y ) that we’ve learned. Actually, if were calculating p ( y | x ) in order to make a prediction, then we don’t actually need to calculate the denominator, since arg max y p ( y | x ) = arg max y p ( x | y ) p ( y ) p ( x ) = arg max y p ( x | y ) p ( y ) . 4.1 Gaussian discriminant analysis The first generative learning algorithm that we’ll look at is Gaussian discrim- inant analysis (GDA). In this model, we’ll assume that p ( x | y ) is distributed according to a multivariate normal distribution. Let’s talk briefly about the properties of multivariate normal distributions before moving on to the GDA model itself. 4.1.1 The multivariate normal distribution The multivariate normal distribution in d -dimensions, also called the multi- variate Gaussian distribution, is parameterized by a mean vector μ ∈ R d and a covariance matrix Σ ∈ R d × d , where Σ ≥ 0 is symmetric and positive semi-definite. Also written “ N ( μ, Σ)”, its density is given by: p ( x ; μ, Σ) = 1 (2 π ) d/ 2 | Σ | 1 / 2 exp - 1 2 ( x - μ ) T Σ - 1 ( x - μ ) . In the equation above, “ | Σ | ” denotes the determinant of the matrix Σ. For a random variable X distributed N ( μ, Σ), the mean is (unsurpris- ingly) given by μ : E[ X ] = Z x x p ( x ; μ, Σ) dx = μ The covariance of a vector-valued random variable Z is defined as Cov( Z ) = E[( Z - E[ Z ])( Z - E[ Z ]) T ]. This generalizes the notion of the variance of a

36 real-valued random variable. The covariance can also be defined as Cov( Z ) = E[ ZZ T ] - (E[ Z ])(E[ Z ]) T . (You should be able to prove to yourself that these two definitions are equivalent.) If X ∼ N ( μ, Σ), then Cov( X ) = Σ . Here are some examples of what the density of a Gaussian distribution looks like: -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 The left-most figure shows a Gaussian with mean zero (that is, the 2x1 zero-vector) and covariance matrix Σ = I (the 2x2 identity matrix). A Gaus- sian with zero mean and identity covariance is also called the standard nor- mal distribution . The middle figure shows the density of a Gaussian with zero mean and Σ = 0 . 6 I ; and in the rightmost figure shows one with , Σ = 2 I . We see that as Σ becomes larger, the Gaussian becomes more “spread-out,” and as it becomes smaller, the distribution becomes more “compressed.” Let’s look at some more examples. -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 The figures above show Gaussians with mean 0, and with covariance matrices respectively Σ = 1 0 0 1 ; Σ = 1 0.5 0.5 1 ; Σ = 1 0.8 0.8 1 . The leftmost figure shows the familiar standard normal distribution, and we see that as we increase the off-diagonal entry in Σ, the density becomes more

37 “compressed” towards the 45 ◦ line (given by x 1 = x 2 ). We can see this more clearly when we look at the contours of the same three densities: -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 Here’s one last set of examples generated by varying Σ: -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 The plots above used, respectively, Σ = 1 -0.5 -0.5 1 ; Σ = 1 -0.8 -0.8 1 ; Σ = 3 0.8 0.8 1 . From the leftmost and middle figures, we see that by decreasing the off- diagonal elements of the covariance matrix, the density now becomes “com- pressed” again, but in the opposite direction. Lastly, as we vary the pa- rameters, more generally the contours will form ellipses (the rightmost figure showing an example). As our last set of examples, fixing Σ = I , by varying μ , we can also move the mean of the density around. -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 0.05 0.1 0.15 0.2 0.25

38 The figures above were generated using Σ = I , and respectively μ = 1 0 ; μ = -0.5 0 ; μ = -1 -1.5 . 4.1.2 The Gaussian discriminant analysis model When we have a classification problem in which the input features x are continuous-valued random variables, we can then use the Gaussian Discrim- inant Analysis (GDA) model, which models p ( x | y ) using a multivariate nor- mal distribution. The model is: y ∼ Bernoulli( φ ) x | y = 0 ∼ N ( μ 0 , Σ) x | y = 1 ∼ N ( μ 1 , Σ) Writing out the distributions, this is: p ( y ) = φ y (1 - φ ) 1 - y p ( x | y = 0) = 1 (2 π ) d/ 2 | Σ | 1 / 2 exp - 1 2 ( x - μ 0 ) T Σ - 1 ( x - μ 0 ) p ( x | y = 1) = 1 (2 π ) d/ 2 | Σ | 1 / 2 exp - 1 2 ( x - μ 1 ) T Σ - 1 ( x - μ 1 ) Here, the parameters of our model are φ , Σ, μ 0 and μ 1 . (Note that while there’re two different mean vectors μ 0 and μ 1 , this model is usually applied using only one covariance matrix Σ.) The log-likelihood of the data is given by ‘ ( φ, μ 0 , μ 1 , Σ) = log n Y i =1 p ( x ( i ) , y ( i ) ; φ, μ 0 , μ 1 , Σ) = log n Y i =1 p ( x ( i ) | y ( i ) ; μ 0 , μ 1 , Σ) p ( y ( i ) ; φ ) .

39 By maximizing ‘ with respect to the parameters, we find the maximum like- lihood estimate of the parameters (see problem set 1) to be: φ = 1 n n X i =1 1 { y ( i ) = 1 } μ 0 = ∑ n i =1 1 { y ( i ) = 0 } x ( i ) ∑ n i =1 1 { y ( i ) = 0 } μ 1 = ∑ n i =1 1 { y ( i ) = 1 } x ( i ) ∑ n i =1 1 { y ( i ) = 1 } Σ = 1 n n X i =1 ( x ( i ) - μ y ( i ) )( x ( i ) - μ y ( i ) ) T . Pictorially, what the algorithm is doing can be seen in as follows: -2 -1 0 1 2 3 4 5 6 7 -7 -6 -5 -4 -3 -2 -1 0 1 Shown in the figure are the training set, as well as the contours of the two Gaussian distributions that have been fit to the data in each of the two classes. Note that the two Gaussians have contours that are the same shape and orientation, since they share a covariance matrix Σ, but they have different means μ 0 and μ 1 . Also shown in the figure is the straight line giving the decision boundary at which p ( y = 1 | x ) = 0 . 5. On one side of the boundary, we’ll predict y = 1 to be the most likely outcome, and on the other side, we’ll predict y = 0.

40 4.1.3 Discussion: GDA and logistic regression The GDA model has an interesting relationship to logistic regression. If we view the quantity p ( y = 1 | x ; φ, μ 0 , μ 1 , Σ) as a function of x , we’ll find that it can be expressed in the form p ( y = 1 | x ; φ, Σ , μ 0 , μ 1 ) = 1 1 + exp( - θ T x ) , where θ is some appropriate function of φ, Σ , μ 0 , μ 1 . 1 This is exactly the form that logistic regression—a discriminative algorithm—used to model p ( y = 1 | x ). When would we prefer one model over another? GDA and logistic regres- sion will, in general, give different decision boundaries when trained on the same dataset. Which is better? We just argued that if p ( x | y ) is multivariate gaussian (with shared Σ), then p ( y | x ) necessarily follows a logistic function. The converse, however, is not true; i.e., p ( y | x ) being a logistic function does not imply p ( x | y ) is multivariate gaussian. This shows that GDA makes stronger modeling as- sumptions about the data than does logistic regression. It turns out that when these modeling assumptions are correct, then GDA will find better fits to the data, and is a better model. Specifically, when p ( x | y ) is indeed gaus- sian (with shared Σ), then GDA is asymptotically efficient . Informally, this means that in the limit of very large training sets (large n ), there is no algorithm that is strictly better than GDA (in terms of, say, how accurately they estimate p ( y | x )). In particular, it can be shown that in this setting, GDA will be a better algorithm than logistic regression; and more generally, even for small training set sizes, we would generally expect GDA to better. In contrast, by making significantly weaker assumptions, logistic regres- sion is also more robust and less sensitive to incorrect modeling assumptions. There are many different sets of assumptions that would lead to p ( y | x ) taking the form of a logistic function. For example, if x | y = 0 ∼ Poisson( λ 0 ), and x | y = 1 ∼ Poisson( λ 1 ), then p ( y | x ) will be logistic. Logistic regression will also work well on Poisson data like this. But if we were to use GDA on such data—and fit Gaussian distributions to such non-Gaussian data—then the results will be less predictable, and GDA may (or may not) do well. To summarize: GDA makes stronger modeling assumptions, and is more data efficient (i.e., requires less training data to learn “well”) when the mod- eling assumptions are correct or at least approximately correct. Logistic 1 This uses the convention of redefining the x ( i ) ’s on the right-hand-side to be ( d + 1)- dimensional vectors by adding the extra coordinate x ( i ) 0 = 1; see problem set 1.

41 regression makes weaker assumptions, and is significantly more robust to deviations from modeling assumptions. Specifically, when the data is in- deed non-Gaussian, then in the limit of large datasets, logistic regression will almost always do better than GDA. For this reason, in practice logistic re- gression is used more often than GDA. (Some related considerations about discriminative vs. generative models also apply for the Naive Bayes algo- rithm that we discuss next, but the Naive Bayes algorithm is still considered a very good, and is certainly also a very popular, classification algorithm.) 4.2 Naive bayes In GDA, the feature vectors x were continuous, real-valued vectors. Let’s now talk about a different learning algorithm in which the x j ’s are discrete- valued. For our motivating example, consider building an email spam filter using machine learning. Here, we wish to classify messages according to whether they are unsolicited commercial (spam) email, or non-spam email. After learning to do this, we can then have our mail reader automatically filter out the spam messages and perhaps place them in a separate mail folder. Classifying emails is one example of a broader set of problems called text classification . Let’s say we have a training set (a set of emails labeled as spam or non- spam). We’ll begin our construction of our spam filter by specifying the features x j used to represent an email. We will represent an email via a feature vector whose length is equal to the number of words in the dictionary. Specifically, if an email contains the j -th word of the dictionary, then we will set x j = 1; otherwise, we let x j = 0. For instance, the vector x =            1 0 0 . . . 1 . . . 0            a aardvark aardwolf . . . buy . . . zygmurgy is used to represent an email that contains the words “a” and “buy,” but not

42 “aardvark,” “aardwolf” or “zygmurgy.” 2 The set of words encoded into the feature vector is called the vocabulary , so the dimension of x is equal to the size of the vocabulary. Having chosen our feature vector, we now want to build a generative model. So, we have to model p ( x | y ). But if we have, say, a vocabulary of 50000 words, then x ∈ { 0 , 1 } 50000 ( x is a 50000-dimensional vector of 0’s and 1’s), and if we were to model x explicitly with a multinomial distribution over the 2 50000 possible outcomes, then we’d end up with a (2 50000 - 1)-dimensional parameter vector. This is clearly too many parameters. To model p ( x | y ), we will therefore make a very strong assumption. We will assume that the x i ’s are conditionally independent given y . This assumption is called the Naive Bayes (NB) assumption , and the resulting algorithm is called the Naive Bayes classifier . For instance, if y = 1 means spam email; “buy” is word 2087 and “price” is word 39831; then we are assuming that if I tell you y = 1 (that a particular piece of email is spam), then knowledge of x 2087 (knowledge of whether “buy” appears in the message) will have no effect on your beliefs about the value of x 39831 (whether “price” appears). More formally, this can be written p ( x 2087 | y ) = p ( x 2087 | y, x 39831 ). (Note that this is not the same as saying that x 2087 and x 39831 are independent, which would have been written “ p ( x 2087 ) = p ( x 2087 | x 39831 )”; rather, we are only assuming that x 2087 and x 39831 are conditionally independent given y .) We now have: p ( x 1 , . . . , x 50000 | y ) = p ( x 1 | y ) p ( x 2 | y, x 1 ) p ( x 3 | y, x 1 , x 2 ) · · · p ( x 50000 | y, x 1 , . . . , x 49999 ) = p ( x 1 | y ) p ( x 2 | y ) p ( x 3 | y ) · · · p ( x 50000 | y ) = d Y j =1 p ( x j | y ) The first equality simply follows from the usual properties of probabilities, and the second equality used the NB assumption. We note that even though 2 Actually, rather than looking through an English dictionary for the list of all English words, in practice it is more common to look through our training set and encode in our feature vector only the words that occur at least once there. Apart from reducing the number of words modeled and hence reducing our computational and space requirements, this also has the advantage of allowing us to model/include as a feature many words that may appear in your email (such as “cs229”) but that you won’t find in a dictionary. Sometimes (as in the homework), we also exclude the very high frequency words (which will be words like “the,” “of,” “and”; these high frequency, “content free” words are called stop words ) since they occur in so many documents and do little to indicate whether an email is spam or non-spam.