Contents
Andrew Ng

Introduction to Supervised Learning

Let's start by talking about a few examples of supervised learning problems. Suppose we have a dataset giving the living areas and prices of 47 47 4747 houses from Portland, Oregon:
Living area ( feet 2 ) feet 2 ("feet"^(2))\left({\text{feet}}^{2}\right) Price ( 1000 $ s ) 1000 $ s (1000$"s")\left(1000\$\text{s}\right)
2104 400
1600 330
2400 369
1416 232
3000 540
⋮ ⋮ vdots\vdots ⋮ ⋮ vdots\vdots
We can plot this data:
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 ) x ( i ) x^((i))x^{(i)} to denote the "input" variables (living area in this example), also called input features, and y ( i ) y ( i ) y^((i))y^{(i)} to denote the "output" or target variable that we are trying to predict (price). A pair ( x ( i ) , y ( i ) ) x ( i ) , y ( i ) (x^((i)),y^((i)))\left(x^{(i)}, y^{(i)}\right) is called a training example, and the dataset that we'll be using to learn–a list of n n nn training examples { ( x ( i ) , y ( i ) ) ; i = 1 , … , n } x ( i ) , y ( i ) ; i = 1 , … , n {(x^((i)),y^((i)));i=1,dots,n}\left\{\left(x^{(i)}, y^{(i)}\right) ; i=1, \ldots, n\right\}–is called a training set. Note that the superscript " i ( i ) i ( i ) i(i)i(i)" in the notation is simply an index into the training set, and has nothing to do with exponentiation. We will also use X X X\mathcal{X} denote the space of input values, and Y Y Y\mathcal{Y} the space of output values. In this example, X = Y = R X = Y = R X=Y=R\mathcal{X}=\mathcal{Y}=\mathbb{R}.
To describe the supervised learning problem slightly more formally, our goal is, given a training set, to learn a function h : X ↦ Y h : X ↦ Y h:X|->Yh: \mathcal{X} \mapsto \mathcal{Y} so that h ( x ) h ( x ) h(x)h(x) is a "good" predictor for the corresponding value of y y yy. For historical reasons, this function h h hh is called a hypothesis. Seen pictorially, the process is therefore like this:
When the target variable that we're trying to predict is continuous, such as in our housing example, we call the learning problem a regression problem. When y y yy can take on only a small number of discrete values (such as if, given the living area, we wanted to predict if a dwelling is a house or an apartment, say), we call it a classification problem.
Part I

Linear Regression

To make our housing example more interesting, let's consider a slightly richer dataset in which we also know the number of bedrooms in each house:
Living area ( feet 2 ) feet 2 ("feet"^(2))\left(\text{feet}^{2}\right) #bedrooms Price ( 1000 $ s ) 1000 $ s (1000$"s")\left(1000\$\text{s}\right)
2104 3 400
1600 3 330
2400 3 369
1416 2 232
3000 4 540
⋮ ⋮ vdots\vdots ⋮ ⋮ vdots\vdots ⋮ ⋮ vdots\vdots
Here, the x x xx's are two-dimensional vectors in R 2 R 2 R^(2)\mathbb{R}^{2}. For instance, x 1 ( i ) x 1 ( i ) x_(1)^((i))x_{1}^{(i)} is the living area of the i i ii-th house in the training set, and x 2 ( i ) x 2 ( i ) x_(2)^((i))x_{2}^{(i)} is its number of bedrooms. (In general, when designing a learning problem, it will be up to you to decide what features to choose, so if you are out in Portland gathering housing data, you might also decide to include other features such as whether each house has a fireplace, the number of bathrooms, and so on. We'll say more about feature selection later, but for now let's take the features as given.)
To perform supervised learning, we must decide how we're going to represent functions/hypotheses h h hh in a computer. As an initial choice, let's say we decide to approximate y y yy as a linear function of x x xx:
h θ ( x ) = θ 0 + θ 1 x 1 + θ 2 x 2 h θ ( x ) = θ 0 + θ 1 x 1 + θ 2 x 2 h_(theta)(x)=theta_(0)+theta_(1)x_(1)+theta_(2)x_(2)h_{\theta}(x)=\theta_{0}+\theta_{1} x_{1}+\theta_{2} x_{2}
Here, the θ i θ i theta_(i)\theta_{i}'s are the parameters (also called weights) parameterizing the space of linear functions mapping from X X X\mathcal{X} to Y Y Y\mathcal{Y}. When there is no risk of confusion, we will drop the θ θ theta\theta subscript in h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x), and write it more simply as h ( x ) h ( x ) h(x)h(x). To simplify our notation, we also introduce the convention of letting x 0 = 1 x 0 = 1 x_(0)=1x_{0}=1 (this is the intercept term), so that
h ( x ) = ∑ i = 0 d θ i x i = θ T x , h ( x ) = ∑ i = 0 d   θ i x i = θ T x , h(x)=sum_(i=0)^(d)theta_(i)x_(i)=theta^(T)x,h(x)=\sum_{i=0}^{d} \theta_{i} x_{i}=\theta^{T} x ,
where on the right-hand side above we are viewing θ θ theta\theta and x x xx both as vectors, and here d d dd is the number of input variables (not counting x 0 x 0 x_(0)x_{0}).
Now, given a training set, how do we pick, or learn, the parameters θ θ theta\theta? One reasonable method seems to be to make h ( x ) h ( x ) h(x)h(x) close to y y yy, at least for the training examples we have. To formalize this, we will define a function that measures, for each value of the θ θ theta\theta's, how close the h ( x ( i ) ) h ( x ( i ) ) h(x^((i)))h(x^{(i)})'s are to the corresponding y ( i ) y ( i ) y^((i))y^{(i)}'s. We define the cost function:
J ( θ ) = 1 2 ∑ i = 1 n ( h θ ( x ( i ) ) − y ( i ) ) 2 . J ( θ ) = 1 2 ∑ i = 1 n   ( h θ ( x ( i ) ) − y ( i ) ) 2 . J(theta)=(1)/(2)sum_(i=1)^(n)(h_(theta)(x^((i)))-y^((i)))^(2).J(\theta)=\frac{1}{2} \sum_{i=1}^{n}(h_{\theta}(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. LMS algorithm

We want to choose θ θ theta\theta so as to minimize J ( θ ) J ( θ ) J(theta)J(\theta). To do so, let's use a search algorithm that starts with some "initial guess" for θ θ theta\theta, and that repeatedly changes θ θ theta\theta to make J ( θ ) J ( θ ) J(theta)J(\theta) smaller, until hopefully we converge to a value of θ θ theta\theta that minimizes J ( θ ) J ( θ ) J(theta)J(\theta). Specifically, let's consider the gradient descent algorithm, which starts with some initial θ θ theta\theta, and repeatedly performs the update:
θ j := θ j − α ∂ ∂ θ j J ( θ ) . θ j := θ j − α ∂ ∂ θ j J ( θ ) . theta_(j):=theta_(j)-alpha(del)/(deltheta_(j))J(theta).\theta_{j}:=\theta_{j}-\alpha \frac{\partial}{\partial \theta_{j}} J(\theta) .
(This update is simultaneously performed for all values of j = 0 , … , d j = 0 , … , d j=0,dots,dj=0, \ldots, d.) Here, α α alpha\alpha is called the learning rate. This is a very natural algorithm that repeatedly takes a step in the direction of steepest decrease of J J JJ.
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 case of if we have only one training example ( x , y ) ( x , y ) (x,y)(x, y), so that we can neglect the sum in the definition of J J JJ. We have:
∂ ∂ θ j J ( θ ) = ∂ ∂ θ j 1 2 ( h θ ( x ) − y ) 2 = 2 ⋅ 1 2 ( h θ ( x ) − y ) ⋅ ∂ ∂ θ j ( h θ ( x ) − y ) = ( h θ ( x ) − y ) ⋅ ∂ ∂ θ j ( ∑ i = 0 d θ i x i − y ) = ( h θ ( x ) − y ) x j ∂ ∂ θ j J ( θ ) = ∂ ∂ θ j 1 2 h θ ( x ) − y 2 = 2 ⋅ 1 2 h θ ( x ) − y ⋅ ∂ ∂ θ j h θ ( x ) − y = h θ ( x ) − y ⋅ ∂ ∂ θ j ∑ i = 0 d   θ i x i − y = h θ ( x ) − y x j {:[(del)/(deltheta_(j))J(theta)=(del)/(deltheta_(j))(1)/(2)(h_(theta)(x)-y)^(2)],[=2*(1)/(2)(h_(theta)(x)-y)*(del)/(deltheta_(j))(h_(theta)(x)-y)],[=(h_(theta)(x)-y)*(del)/(deltheta_(j))(sum_(i=0)^(d)theta_(i)x_(i)-y)],[=(h_(theta)(x)-y)x_(j)]:}\begin{aligned} \frac{\partial}{\partial \theta_{j}} J(\theta) &=\frac{\partial}{\partial \theta_{j}} \frac{1}{2}\left(h_{\theta}(x)-y\right)^{2} \\ &=2 \cdot \frac{1}{2}\left(h_{\theta}(x)-y\right) \cdot \frac{\partial}{\partial \theta_{j}}\left(h_{\theta}(x)-y\right) \\ &=\left(h_{\theta}(x)-y\right) \cdot \frac{\partial}{\partial \theta_{j}}\left(\sum_{i=0}^{d} \theta_{i} x_{i}-y\right) \\ &=\left(h_{\theta}(x)-y\right) x_{j} \end{aligned}
For a single training example, this gives the update rule:[1]
θ j := θ j + α ( y ( i ) − h θ ( x ( i ) ) ) x j ( i ) . θ j := θ j + α y ( i ) − h θ ( x ( i ) ) x j ( i ) . theta_(j):=theta_(j)+alpha(y^((i))-h_(theta)(x^((i))))x_(j)^((i)).\theta_{j}:=\theta_{j}+\alpha\left(y^{(i)}-h_{\theta}(x^{(i)})\right) x_{j}^{(i)} .
The rule is called the LMS update rule (LMS stands for "least mean squares"), and is also known as the Widrow-Hoff learning rule. This rule has several properties that seem natural and intuitive. For instance, the magnitude of the update is proportional to the error term ( y ( i ) − h θ ( x ( i ) ) ) ( y ( i ) − h θ ( x ( i ) ) ) (y^((i))-h_(theta)(x^((i))))(y^{(i)}-h_{\theta}(x^{(i)})); thus, for instance, if we are encountering a training example on which our prediction nearly matches the actual value of y ( i ) y ( i ) y^((i))y^{(i)}, then we find that there is little need to change the parameters; in contrast, a larger change to the parameters will be made if our prediction h θ ( x ( i ) ) h θ ( x ( i ) ) h_(theta)(x^((i)))h_{\theta}(x^{(i)}) has a large error (i.e., if it is very far from y ( i ) y ( i ) y^((i))y^{(i)}).
We'd derived the LMS rule for when there was only a single training example. There are two ways to modify this method for a training set of more than one example. The first is replace it with the following algorithm:
Repeat until convergence { { {\{
(1) θ j := θ j + α ∑ i = 1 n ( y ( i ) − h θ ( x ( i ) ) ) x j ( i ) , ( for every j ) (1) θ j := θ j + α ∑ i = 1 n   y ( i ) − h θ x ( i ) x j ( i ) , ( for every  j ) {:(1)theta_(j):=theta_(j)+alphasum_(i=1)^(n)(y^((i))-h_(theta)(x^((i))))x_(j)^((i))","("for every "j):}\begin{equation} \theta_{j}:=\theta_{j}+\alpha \sum_{i=1}^{n}\left(y^{(i)}-h_{\theta}\left(x^{(i)}\right)\right) x_{j}^{(i)},(\text {for every } j) \end{equation}
} } }\}
By grouping the updates of the coordinates into an update of the vector θ θ theta\theta, we can rewrite update (1) in a slightly more succinct way:
θ := θ + α ∑ i = 1 n ( y ( i ) − h θ ( x ( i ) ) ) x ( i ) θ := θ + α ∑ i = 1 n   y ( i ) − h θ ( x ( i ) ) x ( i ) theta:=theta+alphasum_(i=1)^(n)(y^((i))-h_(theta)(x^((i))))x^((i))\theta:=\theta+\alpha \sum_{i=1}^{n}\left(y^{(i)}-h_{\theta}(x^{(i)})\right) x^{(i)}
The reader can easily verify that the quantity in the summation in the update rule above is just ∂ J ( θ ) / ∂ θ j ∂ J ( θ ) / ∂ θ j del J(theta)//deltheta_(j)\partial J(\theta) / \partial \theta_{j} (for the original definition of J J JJ). So, this is simply gradient descent on the original cost function J J JJ. This method looks at every example in the entire training set on every step, and is called batch gradient descent. Note that, while gradient descent can be susceptible to local minima in general, the optimization problem we have posed here for linear regression has only one global, and no other local, optima; thus gradient descent always converges (assuming the learning rate α α alpha\alpha is not too large) to the global minimum. Indeed, J J JJ is a convex quadratic function. Here is an example of gradient descent as it is run to minimize a quadratic function.
The ellipses shown above are the contours of a quadratic function. Also shown is the trajectory taken by gradient descent, which was initialized at ( 48 , 30 ) ( 48 , 30 ) (48,30)(48,30). The x x xx's in the figure (joined by straight lines) mark the successive values of θ θ theta\theta that gradient descent went through.
When we run batch gradient descent to fit θ θ theta\theta on our previous dataset, to learn to predict housing price as a function of living area, we obtain θ 0 = 71.27 , θ 1 = 0.1345 θ 0 = 71.27 , θ 1 = 0.1345 theta_(0)=71.27,theta_(1)=0.1345\theta_{0}=71.27, \theta_{1}=0.1345. If we plot h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x) as a function of x x xx (area), along with the training data, we obtain the following figure:
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 θ 0 = 89.60 , θ 1 = 0.1392 , θ 2 = − 8.738 theta_(0)=89.60,theta_(1)=0.1392,theta_(2)=-8.738\theta_{0}=89.60, \theta_{1}=0.1392, \theta_{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 { { {\{
quad\quad for i = 1 i = 1 i=1i = 1 to n , { n , { n,{n, \{
(2) θ j := θ j + α ( y ( i ) − h θ ( x ( i ) ) ) x j ( i ) , ( for every j ) (2) θ j := θ j + α y ( i ) − h θ ( x ( i ) ) x j ( i ) , ( for every  j ) {:(2)theta_(j):=theta_(j)+alpha(y^((i))-h_(theta)(x^((i))))x_(j)^((i))","quad("for every "j):}\begin{equation} \theta_{j}:=\theta_{j}+\alpha\left(y^{(i)}-h_{\theta}(x^{(i)})\right)x^{(i)}_{j},\quad(\text{for every }j) \tag{2} \end{equation}
} } quad}\quad \}
} } }\}
By grouping the updates of the coordinates into an update of the vector θ θ theta\theta, we can rewrite update (2) in a slightly more succinct way:
θ := θ + α ( y ( i ) − h θ ( x ( i ) ) ) x ( i ) θ := θ + α y ( i ) − h θ ( x ( i ) ) x ( i ) theta:=theta+alpha(y^((i))-h_(theta)(x^((i))))x^((i))\theta:=\theta+\alpha\left(y^{(i)}-h_{\theta}(x^{(i)})\right) 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 n nn is large–stochastic gradient descent can start making progress right away, and continues to make progress with each example it looks at. Often, stochastic gradient descent gets θ θ theta\theta "close" to the minimum much faster than batch gradient descent. (Note however that it may never "converge" to the minimum, and the parameters θ θ theta\theta will keep oscillating around the minimum of J ( θ ) J ( θ ) J(theta)J(\theta); but in practice most of the values near the minimum will be reasonably good approximations to the true minimum.[2]) For these reasons, particularly when the training set is large, stochastic gradient descent is often preferred over batch gradient descent.

2. The normal equations

Gradient descent gives one way of minimizing J J JJ. Let's discuss a second way of doing so, this time performing the minimization explicitly and without resorting to an iterative algorithm. In this method, we will minimize J J JJ by explicitly taking its derivatives with respect to the θ j θ j theta_(j)\theta_{j}'s, and setting them to zero. To enable us to do this without having to write reams of algebra and pages full of matrices of derivatives, let's introduce some notation for doing calculus with matrices.

2.1. Matrix derivatives

For a function f : R n × d ↦ R f : R n × d ↦ R f:R^(n xx d)|->Rf: \mathbb{R}^{n \times d} \mapsto \mathbb{R} mapping from n n nn-by- d d dd matrices to the real numbers, we define the derivative of f f ff with respect to A A AA to be:
∇ A f ( A ) = [ ∂ f ∂ A 11 ⋯ ∂ f ∂ A 1 d ⋮ ⋱ ⋮ ∂ f ∂ A n 1 ⋯ ∂ f ∂ A n d ] ∇ A f ( A ) = ∂ f ∂ A 11 ⋯ ∂ f ∂ A 1 d ⋮ ⋱ ⋮ ∂ f ∂ A n 1 ⋯ ∂ f ∂ A n d grad_(A)f(A)=[[(del f)/(delA_(11)),cdots,(del f)/(delA_(1d))],[vdots,ddots,vdots],[(del f)/(delA_(n1)),cdots,(del f)/(delA_(nd))]]\nabla_{A} f(A)=\left[\begin{array}{ccc} \frac{\partial f}{\partial A_{11}} & \cdots & \frac{\partial f}{\partial A_{1 d}} \\ \vdots & \ddots & \vdots \\ \frac{\partial f}{\partial A_{n 1}} & \cdots & \frac{\partial f}{\partial A_{n d}} \end{array}\right]
Thus, the gradient ∇ A f ( A ) ∇ A f ( A ) grad_(A)f(A)\nabla_{A} f(A) is itself an n n nn-by- d d dd matrix, whose ( i , j ) ( i , j ) (i,j)(i, j)-element is ∂ f / ∂ A i j ∂ f / ∂ A i j del f//delA_(ij)\partial f / \partial A_{i j}. For example, suppose A = [ A 11 A 12 A 21 A 22 ] A = A 11      A 12 A 21      A 22 A=[[A_(11),A_(12)],[A_(21),A_(22)]]A=\left[\begin{array}{ll}A_{11} & A_{12} \\ A_{21} & A_{22}\end{array}\right] is a 2-by-2 matrix, and the function f : R 2 × 2 ↦ R f : R 2 × 2 ↦ R f:R^(2xx2)|->Rf: \mathbb{R}^{2 \times 2} \mapsto \mathbb{R} is given by
f ( A ) = 3 2 A 11 + 5 A 12 2 + A 21 A 22 . f ( A ) = 3 2 A 11 + 5 A 12 2 + A 21 A 22 . f(A)=(3)/(2)A_(11)+5A_(12)^(2)+A_(21)A_(22).f(A)=\frac{3}{2} A_{11}+5 A_{12}^{2}+A_{21} A_{22} .
Here, A i j A i j A_(ij)A_{i j} denotes the ( i , j ) ( i , j ) (i,j)(i, j) entry of the matrix A A AA. We then have
∇ A f ( A ) = [ 3 2 10 A 12 A 22 A 21 ] . ∇ A f ( A ) = 3 2 10 A 12 A 22 A 21 . grad_(A)f(A)=[[(3)/(2),10A_(12)],[A_(22),A_(21)]].\nabla_{A} f(A)=\left[\begin{array}{cc} \frac{3}{2} & 10 A_{12} \\ A_{22} & A_{21} \end{array}\right] .

2.2. Least squares revisited

Armed with the tools of matrix derivatives, let us now proceed to find in closed-form the value of θ θ theta\theta that minimizes J ( θ ) J ( θ ) J(theta)J(\theta). We begin by re-writing J J JJ in matrix-vectorial notation.
Given a training set, define the design matrix X X XX to be the n n nn-by- d d dd matrix (actually n n nn-by- d + 1 d + 1 d+1d+1, if we include the intercept term) that contains the training examples' input values in its rows:
X = [ − ( x ( 1 ) ) T − − ( x ( 2 ) ) T − ⋮ − ( x ( n ) ) T − ] . X = − x ( 1 ) T − − x ( 2 ) T − ⋮ − x ( n ) T − . X=[[-(x^((1)))^(T)-],[-(x^((2)))^(T)-],[vdots],[-(x^((n)))^(T)-]].X=\left[\begin{array}{c} -\left(x^{(1)}\right)^{T}- \\ -\left(x^{(2)}\right)^{T}- \\ \vdots \\ -\left(x^{(n)}\right)^{T}- \end{array}\right].
Also, let y → y → vec(y)\vec{y} be the n n nn-dimensional vector containing all the target values from the training set:
y → = [ y ( 1 ) y ( 2 ) ⋮ y ( n ) ] . y → = y ( 1 ) y ( 2 ) ⋮ y ( n ) . vec(y)=[[y^((1))],[y^((2))],[vdots],[y^((n))]].\vec{y}=\left[\begin{array}{c} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(n)} \end{array}\right].
Now, since h θ ( x ( i ) ) = ( x ( i ) ) T θ h θ x ( i ) = x ( i ) T θ h_(theta)(x^((i)))=(x^((i)))^(T)thetah_{\theta}\left(x^{(i)}\right)=\left(x^{(i)}\right)^{T} \theta, we can easily verify that
X θ − y → = [ ( x ( 1 ) ) T θ ⋮ ( x ( n ) ) T θ ] − [ y ( 1 ) ⋮ y ( n ) ] = [ h θ ( x ( 1 ) ) − y ( 1 ) ⋮ h θ ( x ( n ) ) − y ( n ) ] . X θ − y → = x ( 1 ) T θ ⋮ x ( n ) T θ − y ( 1 ) ⋮ y ( n ) = h θ x ( 1 ) − y ( 1 ) ⋮ h θ x ( n ) − y ( n ) . {:[X theta- vec(y)=[[(x^((1)))^(T)theta],[vdots],[(x^((n)))^(T)theta]]-[[y^((1))],[vdots],[y^((n))]]],[=[[h_(theta)(x^((1)))-y^((1))],[vdots],[h_(theta)(x^((n)))-y^((n))]].]:}\begin{aligned} X \theta-\vec{y} &=\left[\begin{array}{c} \left(x^{(1)}\right)^{T} \theta \\ \vdots \\ \left(x^{(n)}\right)^{T} \theta \end{array}\right]-\left[\begin{array}{c} y^{(1)} \\ \vdots \\ y^{(n)} \end{array}\right] \\ &=\left[\begin{array}{c} h_{\theta}\left(x^{(1)}\right)-y^{(1)} \\ \vdots \\ h_{\theta}\left(x^{(n)}\right)-y^{(n)} \end{array}\right]. \end{aligned}
Thus, using the fact that for a vector z z zz, we have that z T z = ∑ i z i 2 z T z = ∑ i   z i 2 z^(T)z=sum_(i)z_(i)^(2)z^{T} z=\sum_{i} z_{i}^{2}:
1 2 ( X θ − y → ) T ( X θ − y → ) = 1 2 ∑ i = 1 n ( h θ ( x ( i ) ) − y ( i ) ) 2 = J ( θ ) 1 2 ( X θ − y → ) T ( X θ − y → ) = 1 2 ∑ i = 1 n   ( h θ ( x ( i ) ) − y ( i ) ) 2 = J ( θ ) {:[(1)/(2)(X theta- vec(y))^(T)(X theta- vec(y))=(1)/(2)sum_(i=1)^(n)(h_(theta)(x^((i)))-y^((i)))^(2)],[=J(theta)]:}\begin{aligned} \frac{1}{2}(X \theta-\vec{y})^{T}(X \theta-\vec{y}) &=\frac{1}{2} \sum_{i=1}^{n}(h_{\theta}(x^{(i)})-y^{(i)})^{2} \\ &=J(\theta) \end{aligned}
Finally, to minimize J J JJ, let's find its derivatives with respect to θ θ theta\theta. 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 → ∇ θ 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 → {:[grad_(theta)J(theta)=grad_(theta)(1)/(2)(X theta- vec(y))^(T)(X theta- vec(y))],[=(1)/(2)grad_(theta)((X theta)^(T)X theta-(X theta)^(T)( vec(y))- vec(y)^(T)(X theta)+ vec(y)^(T)( vec(y)))],[=(1)/(2)grad_(theta)(theta^(T)(X^(T)X)theta- vec(y)^(T)(X theta)- vec(y)^(T)(X theta))],[=(1)/(2)grad_(theta)(theta^(T)(X^(T)X)theta-2(X^(T)( vec(y)))^(T)theta)],[=(1)/(2)(2X^(T)X theta-2X^(T)( vec(y)))],[=X^(T)X theta-X^(T) vec(y)]:}\begin{aligned} \nabla_{\theta} J(\theta) &=\nabla_{\theta} \frac{1}{2}(X \theta-\vec{y})^{T}(X \theta-\vec{y}) \\ &=\frac{1}{2} \nabla_{\theta}\left((X \theta)^{T} X \theta-(X \theta)^{T} \vec{y}-\vec{y}^{T}(X \theta)+\vec{y}^{T} \vec{y}\right) \\ &=\frac{1}{2} \nabla_{\theta}\left(\theta^{T}\left(X^{T} X\right) \theta-\vec{y}^{T}(X \theta)-\vec{y}^{T}(X \theta)\right) \\ &=\frac{1}{2} \nabla_{\theta}\left(\theta^{T}\left(X^{T} X\right) \theta-2\left(X^{T} \vec{y}\right)^{T} \theta\right) \\ &=\frac{1}{2}\left(2 X^{T} X \theta-2 X^{T} \vec{y}\right) \\ &=X^{T} X \theta-X^{T} \vec{y} \end{aligned}
In the third step, we used the fact that a T b = b T a a T b = b T a a^(T)b=b^(T)aa^{T} b=b^{T} a, and in the fifth step used the facts ∇ x b T x = b ∇ x b T x = b grad_(x)b^(T)x=b\nabla_{x} b^{T} x=b and ∇ x x T A x = 2 A x ∇ x x T A x = 2 A x grad_(x)x^(T)Ax=2Ax\nabla_{x} x^{T} A x=2 A x for symmetric matrix A A AA (for more details, see Section 4.3 of "Linear Algebra Review and Reference"). To minimize J J JJ, we set its derivatives to zero, and obtain the normal equations:
X T X θ = X T y → X T X θ = X T y → X^(T)X theta=X^(T) vec(y)X^{T} X \theta=X^{T} \vec{y}
Thus, the value of θ θ theta\theta that minimizes J ( θ ) J ( θ ) J(theta)J(\theta) is given in closed form by the equation
θ = ( X T X ) − 1 X T y → . θ = X T X − 1 X T y → . theta=(X^(T)X)^(-1)X^(T) vec(y).\begin{equation*} \theta=\left(X^{T} X\right)^{-1} X^{T} \vec{y}. \end{equation*}[3]

3. Probabilistic interpretation

When faced with a regression problem, why might linear regression, and specifically why might the least-squares cost function J J JJ, 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 ) , y ( i ) = θ T x ( i ) + ϵ ( i ) , y^((i))=theta^(T)x^((i))+epsilon^((i)),y^{(i)}=\theta^{T} x^{(i)}+\epsilon^{(i)},
where ϵ ( i ) ϵ ( i ) epsilon^((i))\epsilon^{(i)} is an error term that captures either unmodeled effects (such as if there are some features very pertinent to predicting housing price, but that we'd left out of the regression), or random noise. Let us further assume that the ϵ ( i ) ϵ ( i ) epsilon^((i))\epsilon^{(i)} are distributed IID (independently and identically distributed) according to a Gaussian distribution (also called a Normal distribution) with mean zero and some variance σ 2 σ 2 sigma^(2)\sigma^{2}. We can write this assumption as " ϵ ( i ) ∼ ϵ ( i ) ∼ epsilon^((i))∼\epsilon^{(i)} \sim N ( 0 , σ 2 ) N 0 , σ 2 N(0,sigma^(2))\mathcal{N}\left(0, \sigma^{2}\right)." I.e., the density of ϵ ( i ) ϵ ( i ) epsilon^((i))\epsilon^{(i)} is given by
p ( ϵ ( i ) ) = 1 2 π σ exp ⁡ ( − ( ϵ ( i ) ) 2 2 σ 2 ) . p ( ϵ ( i ) ) = 1 2 π σ exp ⁡ − ϵ ( i ) 2 2 σ 2 . p(epsilon^((i)))=(1)/(sqrt(2pi)sigma)exp(-((epsilon^((i)))^(2))/(2sigma^(2))).p(\epsilon^{(i)})=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(\epsilon^{(i)}\right)^{2}}{2 \sigma^{2}}\right).
This implies that
p ( y ( i ) | x ( i ) ; θ ) = 1 2 π σ exp ⁡ ( − ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) . p ( y ( i ) | x ( i ) ; θ ) = 1 2 π σ exp ⁡ − y ( i ) − θ T x ( i ) 2 2 σ 2 . p(y^((i))|x^((i));theta)=(1)/(sqrt(2pi)sigma)exp(-((y^((i))-theta^(T)x^((i)))^(2))/(2sigma^(2))).p(y^{(i)} | x^{(i)} ; \theta)=\frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}\right) .
The notation " p ( y ( i ) | x ( i ) ; θ ) p y ( i ) | x ( i ) ; θ p(y^((i))|x^((i));theta)p\left(y^{(i)} | x^{(i)} ; \theta\right)" indicates that this is the distribution of y ( i ) y ( i ) y^((i))y^{(i)} given x ( i ) x ( i ) x^((i))x^{(i)} and parameterized by θ θ theta\theta. Note that we should not condition on θ θ theta\theta (" p ( y ( i ) | x ( i ) , θ ) p y ( i ) | x ( i ) , θ p(y^((i))|x^((i)),theta)p\left(y^{(i)} | x^{(i)}, \theta\right)"), since θ θ theta\theta is not a random variable. We can also write the distribution of y ( i ) y ( i ) y^((i))y^{(i)} as y ( i ) | x ( i ) ; θ ∼ N ( θ T x ( i ) , σ 2 ) y ( i ) | x ( i ) ; θ ∼ N θ T x ( i ) , σ 2 y^((i))|x^((i));theta∼N(theta^(T)x^((i)),sigma^(2))y^{(i)} | x^{(i)} ; \theta \sim \mathcal{N}\left(\theta^{T} x^{(i)}, \sigma^{2}\right).
Given X X XX (the design matrix, which contains all the x ( i ) x ( i ) x^((i))x^{(i)}'s) and θ θ theta\theta, what is the distribution of the y ( i ) y ( i ) y^((i))y^{(i)}'s? The probability of the data is given by p ( y → | X ; θ ) p ( y → | X ; θ ) p( vec(y)|X;theta)p(\vec{y} | X ; \theta). This quantity is typically viewed a function of y → y → vec(y)\vec{y} (and perhaps X X XX), for a fixed value of θ θ theta\theta. When we wish to explicitly view this as a function of θ θ theta\theta, we will instead call it the likelihood function:
L ( θ ) = L ( θ ; X , y → ) = p ( y → ∣ X ; θ ) . L ( θ ) = L ( θ ; X , y → ) = p ( y → ∣ X ; θ ) . L(theta)=L(theta;X, vec(y))=p( vec(y)∣X;theta).L(\theta)=L(\theta ; X, \vec{y})=p(\vec{y} \mid X ; \theta).
Note that by the independence assumption on the ϵ ( i ) ϵ ( i ) epsilon^((i))\epsilon^{(i)}'s (and hence also the y ( i ) y ( i ) y^((i))y^{(i)}'s given the x ( i ) x ( i ) x^((i))x^{(i)}'s), this can also be written
L ( θ ) = ∏ i = 1 n p ( y ( i ) ∣ x ( i ) ; θ ) = ∏ i = 1 n 1 2 π σ exp ⁡ ( − ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) . L ( θ ) = ∏ i = 1 n   p y ( i ) ∣ x ( i ) ; θ = ∏ i = 1 n   1 2 π σ exp ⁡ − y ( i ) − θ T x ( i ) 2 2 σ 2 . {:[L(theta)=prod_(i=1)^(n)p(y^((i))∣x^((i));theta)],[=prod_(i=1)^(n)(1)/(sqrt(2pi)sigma)exp(-((y^((i))-theta^(T)x^((i)))^(2))/(2sigma^(2))).]:}\begin{aligned} L(\theta) &=\prod_{i=1}^{n} p\left(y^{(i)} \mid x^{(i)} ; \theta\right) \\ &=\prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}\right) . \end{aligned}
Now, given this probabilistic model relating the y ( i ) y ( i ) y^((i))y^{(i)}'s and the x ( i ) x ( i ) x^((i))x^{(i)}'s, what is a reasonable way of choosing our best guess of the parameters θ θ theta\theta? The principal of maximum likelihood says that we should choose θ θ theta\theta so as to make the data as high probability as possible. I.e., we should choose θ θ theta\theta to maximize L ( θ ) L ( θ ) L(theta)L(\theta).
Instead of maximizing L ( θ ) L ( θ ) L(theta)L(\theta), we can also maximize any strictly increasing function of L ( θ ) L ( θ ) L(theta)L(\theta). In particular, the derivations will be a bit simpler if we instead maximize the log likelihood ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta):
ℓ ( θ ) = log ⁡ L ( θ ) = log ⁡ ∏ i = 1 n 1 2 π σ exp ⁡ ( − ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) = ∑ i = 1 n log ⁡ 1 2 π σ exp ⁡ ( − ( y ( i ) − θ T x ( i ) ) 2 2 σ 2 ) = n log ⁡ 1 2 π σ − 1 σ 2 ⋅ 1 2 ∑ i = 1 n ( y ( i ) − θ T x ( i ) ) 2 . ℓ ( θ ) = log ⁡ L ( θ ) = log ⁡ ∏ i = 1 n   1 2 π σ exp ⁡ − y ( i ) − θ T x ( i ) 2 2 σ 2 = ∑ i = 1 n   log ⁡ 1 2 π σ exp ⁡ − y ( i ) − θ T x ( i ) 2 2 σ 2 = n log ⁡ 1 2 π σ − 1 σ 2 ⋅ 1 2 ∑ i = 1 n   ( y ( i ) − θ T x ( i ) ) 2 . {:[ℓ(theta)=log L(theta)],[=log prod_(i=1)^(n)(1)/(sqrt(2pi)sigma)exp(-((y^((i))-theta^(T)x^((i)))^(2))/(2sigma^(2)))],[=sum_(i=1)^(n)log (1)/(sqrt(2pi)sigma)exp(-((y^((i))-theta^(T)x^((i)))^(2))/(2sigma^(2)))],[=n log (1)/(sqrt(2pi)sigma)-(1)/(sigma^(2))*(1)/(2)sum_(i=1)^(n)(y^((i))-theta^(T)x^((i)))^(2).]:}\begin{aligned} \ell(\theta) &=\log L(\theta) \\ &=\log \prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}\right) \\ &=\sum_{i=1}^{n} \log \frac{1}{\sqrt{2 \pi} \sigma} \exp \left(-\frac{\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2}}{2 \sigma^{2}}\right) \\ &=n \log \frac{1}{\sqrt{2 \pi} \sigma}-\frac{1}{\sigma^{2}} \cdot \frac{1}{2} \sum_{i=1}^{n}(y^{(i)}-\theta^{T} x^{(i)})^{2} . \end{aligned}
Hence, maximizing ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta) gives the same answer as minimizing
1 2 ∑ i = 1 n ( y ( i ) − θ T x ( i ) ) 2 , 1 2 ∑ i = 1 n   y ( i ) − θ T x ( i ) 2 , (1)/(2)sum_(i=1)^(n)(y^((i))-theta^(T)x^((i)))^(2),\frac{1}{2} \sum_{i=1}^{n}\left(y^{(i)}-\theta^{T} x^{(i)}\right)^{2},
which we recognize to be J ( θ ) J ( θ ) J(theta)J(\theta), our original least-squares cost function.
To summarize: Under the previous probabilistic assumptions on the data, least-squares regression corresponds to finding the maximum likelihood estimate of θ θ theta\theta. This is thus one set of assumptions under which least-squares regression can be justified as a very natural method that's just doing maximum likelihood estimation. (Note however that the probabilistic assumptions are by no means necessary for least-squares to be a perfectly good and rational procedure, and there may–and indeed there are–other natural assumptions that can also be used to justify it.)
Note also that, in our previous discussion, our final choice of θ θ theta\theta did not depend on what was σ 2 σ 2 sigma^(2)\sigma^{2}, and indeed we'd have arrived at the same result even if σ 2 σ 2 sigma^(2)\sigma^{2} were unknown. We will use this fact again later, when we talk about the exponential family and generalized linear models.

4. Locally weighted linear regression

Consider the problem of predicting y y yy from x ∈ R x ∈ R x inRx \in \mathbb{R}. The leftmost figure below shows the result of fitting a y = θ 0 + θ 1 x y = θ 0 + θ 1 x y=theta_(0)+theta_(1)xy=\theta_{0}+\theta_{1} x to a dataset. We see that the data doesn't really lie on straight line, and so the fit is not very good.
Instead, if we had added an extra feature x 2 x 2 x^(2)x^{2}, and fit y = θ 0 + θ 1 x + θ 2 x 2 y = θ 0 + θ 1 x + θ 2 x 2 y=theta_(0)+theta_(1)x+theta_(2)x^(2)y=\theta_{0}+\theta_{1} x+\theta_{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 = ∑ j = 0 5 θ j x j y = ∑ j = 0 5   θ j x j y=sum_(j=0)^(5)theta_(j)x^(j)y=\sum_{j=0}^{5} \theta_{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 ) ( y ) (y)(y) for different living areas ( x ) ( x ) (x)(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 automatically choosing a good set of features.) In this section, let us briefly talk about the locally weighted linear regression (LWR) algorithm which, assuming 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 x xx (i.e., to evaluate h ( x ) h ( x ) h(x)h(x)), we would:
  1. Fit θ θ theta\theta to minimize ∑ i ( y ( i ) − θ T x ( i ) ) 2 ∑ i   ( y ( i ) − θ T x ( i ) ) 2 sum_(i)(y^((i))-theta^(T)x^((i)))^(2)\sum_{i}(y^{(i)}-\theta^{T}x^{(i)})^{2}.
  2. Output θ T x θ T x theta^(T)x\theta^{T}x.
In contrast, the locally weighted linear regression algorithm does the following:
  1. Fit θ θ theta\theta to minimize ∑ i w ( i ) ( y ( i ) − θ T x ( i ) ) 2 ∑ i   w ( i ) ( y ( i ) − θ T x ( i ) ) 2 sum_(i)w^((i))(y^((i))-theta^(T)x^((i)))^(2)\sum_{i}w^{(i)}(y^{(i)}-\theta^{T}x^{(i)})^{2}.
  2. Output θ T x θ T x theta^(T)x\theta^{T}x.
Here, the w ( i ) w ( i ) w^((i))w^{(i)}'s are non-negative valued weights. Intuitively, if w ( i ) w ( i ) w^((i))w^{(i)} is large for a particular value of i i ii, then in picking θ θ theta\theta, we'll try hard to make ( y ( i ) − θ T x ( i ) ) 2 ( y ( i ) − θ T x ( i ) ) 2 (y^((i))-theta^(T)x^((i)))^(2)(y^{(i)}-\theta^{T}x^{(i)})^{2} small. If w ( i ) w ( i ) w^((i))w^{(i)} is small, then the ( y ( i ) − θ T x ( i ) ) 2 ( y ( i ) − θ T x ( i ) ) 2 (y^((i))-theta^(T)x^((i)))^(2)(y^{(i)}-\theta^{T}x^{(i)})^{2} error term will be pretty much ignored in the fit.
A fairly standard choice for the weights is[4]
w ( i ) = exp ⁡ ( − ( x ( i ) − x ) 2 2 τ 2 ) w ( i ) = exp ⁡ − ( x ( i ) − x ) 2 2 τ 2 w^((i))=exp(-((x^((i))-x)^(2))/(2tau^(2)))w^{(i)}=\exp\left(-\frac{(x^{(i)}-x)^{2}}{2\tau^{2}}\right)
Note that the weights depend on the particular point x x xx at which we're trying to evaluate x x xx. Moreover, if | x ( i ) − x | | x ( i ) − x | |x^((i))-x||x^{(i)}-x| is small, then w ( i ) w ( i ) w^((i))w^{(i)} is close to 1 1 11; and if | x ( i ) − x | | x ( i ) − x | |x^((i))-x||x^{(i)}-x| is large, then w ( i ) w ( i ) w^((i))w^{(i)} is small. Hence, θ θ theta\theta is chosen giving a much higher "weight" to the (errors on) training examples close to the query point x x xx. (Note also that while the formula for the weights takes a form that is cosmetically similar to the density of a Gaussian distribution, the w ( i ) w ( i ) w^((i))w^{(i)}'s do not directly have anything to do with Gaussians, and in particular the w ( i ) w ( i ) w^((i))w^{(i)} are not random variables, normally distributed or otherwise.) The parameter τ τ tau\tau controls how quickly the weight of a training example falls off with distance of its x ( i ) x ( i ) x^((i))x^{(i)} from the query point x ; τ x ; τ x;taux ; \tau is called the bandwidth parameter, and is also something that you'll get to experiment with in your homework.
Locally weighted linear regression is the first example we're seeing of a non-parametric algorithm. The (unweighted) linear regression algorithm that we saw earlier is known as a parametric learning algorithm, because it has a fixed, finite number of parameters (the θ i θ i theta_(i)\theta_{i}'s), which are fit to the data. Once we've fit the θ i θ i theta_(i)\theta_{i}'s and stored them away, we no longer need to keep the training data around to make future predictions. In contrast, to make predictions using locally weighted linear regression, we need to keep the entire training set around. The term "non-parametric" (roughly) refers to the fact that the amount of stuff we need to keep in order to represent the hypothesis h h hh grows linearly with the size of the training set.
Part II

Classification and logistic regression

Let's now talk about the classification problem. This is just like the regression problem, except that the values y y yy we now want to predict take on only a small number of discrete values. For now, we will focus on the binary classification problem in which y y yy can take on only two values, 0 0 00 and 1 1 11. (Most of what we say here will also generalize to the multiple-class case.) For instance, if we are trying to build a spam classifier for email, then x ( i ) x ( i ) x^((i))x^{(i)} may be some features of a piece of email, and y y yy may be 1 1 11 if it is a piece of spam mail, and 0 0 00 otherwise. 0 0 00 is also called the negative class, and 1 1 11 the positive class, and they are sometimes also denoted by the symbols "-" and "+." Given x ( i ) x ( i ) x^((i))x^{(i)}, the corresponding y ( i ) y ( i ) y^((i))y^{(i)} is also called the label for the training example.

5. Logistic regression

We could approach the classification problem ignoring the fact that y y yy is discrete-valued, and use our old linear regression algorithm to try to predict y y yy given x x xx. However, it is easy to construct examples where this method performs very poorly. Intuitively, it also doesn't make sense for h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x) to take values larger than 1 1 11 or smaller than 0 0 00 when we know that y ∈ { 0 , 1 } y ∈ { 0 , 1 } y in{0,1}y \in\{0,1\}.
To fix this, let's change the form for our hypotheses h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x). We will choose
h θ ( x ) = g ( θ T x ) = 1 1 + e − θ T x , h θ ( x ) = g θ T x = 1 1 + e − θ T x , h_(theta)(x)=g(theta^(T)x)=(1)/(1+e^(-theta^(T)x)),h_{\theta}(x)=g\left(\theta^{T} x\right)=\frac{1}{1+e^{-\theta^{T} x}},
where
g ( z ) = 1 1 + e − z g ( z ) = 1 1 + e − z g(z)=(1)/(1+e^(-z))g(z)=\frac{1}{1+e^{-z}}
is called the logistic function or the sigmoid function. Here is a plot showing g ( z ) g ( z ) g(z)g(z):
Notice that g ( z ) g ( z ) g(z)g(z) tends towards 1 1 11 as z → ∞ z → ∞ z rarr ooz \rightarrow \infty, and g ( z ) g ( z ) g(z)g(z) tends towards 0 0 00 as z → − ∞ z → − ∞ z rarr-ooz \rightarrow-\infty. Moreover, g ( z ) g ( z ) g(z)\mathrm{g}(\mathrm{z}), and hence also h ( x ) h ( x ) h(x)h(x), is always bounded between 0 0 00 and 1 1 11. As before, we are keeping the convention of letting x 0 = 1 x 0 = 1 x_(0)=1x_{0}=1, so that θ T x = θ 0 + ∑ j = 1 d θ j x j . θ T x = θ 0 + ∑ j = 1 d   θ j x j . theta^(T)x=theta_(0)+sum_(j=1)^(d)theta_(j)x_(j).\theta^{T} x=\theta_{0}+\sum_{j=1}^{d} \theta_{j} x_{j} .
For now, let's take the choice of g g gg as given. Other functions that smoothly increase from 0 0 00 to 1 1 11 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 ′ g ′ g^(')g^{\prime}:
g ′ ( z ) = d d z 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 ) ) . g ′ ( z ) = d d z 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 ) ) . {:[g^(')(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)).]:}\begin{aligned} g^{\prime}(z)&=\frac{d}{d z} \frac{1}{1+e^{-z}} \\ &=\frac{1}{\left(1+e^{-z}\right)^{2}}\left(e^{-z}\right) \\ &=\frac{1}{\left(1+e^{-z}\right)} \cdot\left(1-\frac{1}{\left(1+e^{-z}\right)}\right) \\ &=g(z)(1-g(z)) . \end{aligned}
So, given the logistic regression model, how do we fit θ θ theta\theta for it? Following how we saw least squares regression could be derived as the maximum likelihood 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.
Let us assume that
P ( y = 1 ∣ x ; θ ) = h θ ( x ) P ( y = 0 ∣ x ; θ ) = 1 − h θ ( x ) P ( y = 1 ∣ x ; θ ) = h θ ( x ) P ( y = 0 ∣ x ; θ ) = 1 − h θ ( x ) {:[P(y=1∣x;theta)=h_(theta)(x)],[P(y=0∣x;theta)=1-h_(theta)(x)]:}\begin{aligned} &P(y=1 \mid x ; \theta)=h_{\theta}(x) \\ &P(y=0 \mid x ; \theta)=1-h_{\theta}(x) \end{aligned}
Note that this can be written more compactly as
p ( y ∣ x ; θ ) = ( h θ ( x ) ) y ( 1 − h θ ( x ) ) 1 − y p ( y ∣ x ; θ ) = h θ ( x ) y 1 − h θ ( x ) 1 − y p(y∣x;theta)=(h_(theta)(x))^(y)(1-h_(theta)(x))^(1-y)p(y \mid x ; \theta)=\left(h_{\theta}(x)\right)^{y}\left(1-h_{\theta}(x)\right)^{1-y}
Assuming that the n n nn training examples were generated independently, we can then write down the likelihood of the parameters as
L ( θ ) = p ( y → ∣ X ; θ ) = ∏ i = 1 n p ( y ( i ) ∣ x ( i ) ; θ ) = ∏ i = 1 n ( h θ ( x ( i ) ) ) y ( i ) ( 1 − h θ ( x ( i ) ) ) 1 − y ( i ) L ( θ ) = p ( y → ∣ X ; θ ) = ∏ i = 1 n   p ( y ( i ) ∣ x ( i ) ; θ ) = ∏ i = 1 n   h θ ( x ( i ) ) y ( i ) 1 − h θ ( x ( i ) ) 1 − y ( i ) {:[L(theta)=p( vec(y)∣X;theta)],[=prod_(i=1)^(n)p(y^((i))∣x^((i));theta)],[=prod_(i=1)^(n)(h_(theta)(x^((i))))^(y^((i)))(1-h_(theta)(x^((i))))^(1-y^((i)))]:}\begin{aligned} L(\theta) &=p(\vec{y} \mid X ; \theta) \\ &=\prod_{i=1}^{n} p(y^{(i)} \mid x^{(i)} ; \theta) \\ &=\prod_{i=1}^{n}\left(h_{\theta}(x^{(i)})\right)^{y^{(i)}}\left(1-h_{\theta}(x^{(i)})\right)^{1-y^{(i)}} \end{aligned}
As before, it will be easier to maximize the log likelihood:
ℓ ( θ ) = log ⁡ L ( θ ) = ∑ i = 1 n y ( i ) log ⁡ h ( x ( i ) ) + ( 1 − y ( i ) ) log ⁡ ( 1 − h ( x ( i ) ) ) ℓ ( θ ) = log ⁡ L ( θ ) = ∑ i = 1 n   y ( i ) log ⁡ h ( x ( i ) ) + ( 1 − y ( i ) ) log ⁡ ( 1 − h ( x ( i ) ) ) {:[ℓ(theta)=log L(theta)],[=sum_(i=1)^(n)y^((i))log h(x^((i)))+(1-y^((i)))log(1-h(x^((i))))]:}\begin{aligned} \ell(\theta) &=\log L(\theta) \\ &=\sum_{i=1}^{n} y^{(i)} \log h(x^{(i)})+(1-y^{(i)}) \log (1-h(x^{(i)})) \end{aligned}
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 θ := θ + α ∇ θ ℓ ( θ ) θ := θ + α ∇ θ ℓ ( θ ) theta:=theta+alphagrad_(theta)ℓ(theta)\theta:=\theta+\alpha \nabla_{\theta} \ell(\theta). (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 ) ( x , y ) (x,y)(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 ∂ ∂ θ 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 {:[(del)/(deltheta_(j))ℓ(theta)=(y(1)/(g(theta^(T)x))-(1-y)(1)/(1-g(theta^(T)x)))(del)/(deltheta_(j))g(theta^(T)x)],[=(y(1)/(g(theta^(T)x))-(1-y)(1)/(1-g(theta^(T)x)))g(theta^(T)x)(1-g(theta^(T)x))(del)/(deltheta_(j))theta^(T)x],[=(y(1-g(theta^(T)x))-(1-y)g(theta^(T)x))x_(j)],[=(y-h_(theta)(x))x_(j)]:}\begin{aligned} \frac{\partial}{\partial \theta_{j}} \ell(\theta) &=\left(y \frac{1}{g\left(\theta^{T} x\right)}-(1-y) \frac{1}{1-g\left(\theta^{T} x\right)}\right) \frac{\partial}{\partial \theta_{j}} g\left(\theta^{T} x\right) \\ &=\left(y \frac{1}{g\left(\theta^{T} x\right)}-(1-y) \frac{1}{1-g\left(\theta^{T} x\right)}\right) g\left(\theta^{T} x\right)\left(1-g\left(\theta^{T} x\right)\right) \frac{\partial}{\partial \theta_{j}} \theta^{T} x \\ &=\left(y\left(1-g\left(\theta^{T} x\right)\right)-(1-y) g\left(\theta^{T} x\right)\right) x_{j} \\ &=\left(y-h_{\theta}(x)\right) x_{j} \end{aligned}
Above, we used the fact that g ′ ( z ) = g ( z ) ( 1 − g ( z ) ) g ′ ( z ) = g ( z ) ( 1 − g ( z ) ) g^(')(z)=g(z)(1-g(z))g^{\prime}(z)=g(z)(1-g(z)). This therefore gives us the stochastic gradient ascent rule
θ j := θ j + α ( y ( i ) − h θ ( x ( i ) ) ) x j ( i ) θ j := θ j + α y ( i ) − h θ x ( i ) x j ( i ) theta_(j):=theta_(j)+alpha(y^((i))-h_(theta)(x^((i))))x_(j)^((i))\theta_{j}:=\theta_{j}+\alpha\left(y^{(i)}-h_{\theta}\left(x^{(i)}\right)\right) x_{j}^{(i)}
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 ) ) h θ x ( i ) h_(theta)(x^((i)))h_{\theta}\left(x^{(i)}\right) is now defined as a non-linear function of θ T x ( i ) θ T x ( i ) theta^(T)x^((i))\theta^{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.

6. Digression: The perceptron learning algorithm

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 0 00 or 1 1 11 exactly. To do so, it seems natural to change the definition of g g gg to be the threshold function:
g ( z ) = { 1 if z ≥ 0 0 if z < 0 g ( z ) = 1      if  z ≥ 0 0      if  z < 0 g(z)={[1,"if "z >= 0],[0,"if "z < 0]:}g(z)= \begin{cases}1 & \text {if } z \geq 0 \\ 0 & \text {if } z<0\end{cases}
If we then let h θ ( x ) = g ( θ T x ) h θ ( x ) = g θ T x h_(theta)(x)=g(theta^(T)x)h_{\theta}(x)=g\left(\theta^{T} x\right) as before but using this modified definition of g g gg, and if we use the update rule
θ j := θ j + α ( y ( i ) − h θ ( x ( i ) ) ) x j ( i ) . θ j := θ j + α y ( i ) − h θ x ( i ) x j ( i ) . theta_(j):=theta_(j)+alpha(y^((i))-h_(theta)(x^((i))))x_(j)^((i)).\theta_{j}:=\theta_{j}+\alpha\left(y^{(i)}-h_{\theta}\left(x^{(i)}\right)\right) x_{j}^{(i)} .
then we have the perceptron learning algorithm.
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 predictions with meaningful probabilistic interpretations, or derive the perceptron as a maximum likelihood estimation algorithm.

7. Another algorithm for maximizing ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta)

Returning to logistic regression with g ( z ) g ( z ) g(z)g(z) being the sigmoid function, let's now talk about a different algorithm for maximizing ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta).
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 ↦ R f : R ↦ R f:R|->Rf: \mathbb{R} \mapsto \mathbb{R}, and we wish to find a value of θ θ theta\theta so that f ( θ ) = 0 f ( θ ) = 0 f(theta)=0f(\theta)=0. Here, θ ∈ R θ ∈ R theta inR\theta \in \mathbb{R} is a real number. Newton's method performs the following update:
θ := θ − f ( θ ) f ′ ( θ ) . θ := θ − f ( θ ) f ′ ( θ ) . theta:=theta-(f(theta))/(f^(')(theta)).\theta:=\theta-\frac{f(\theta)}{f^{\prime}(\theta)}.
This method has a natural interpretation in which we can think of it as approximating the function f f ff via a linear function that is tangent to f f ff at the current guess θ θ theta\theta, solving for where that linear function equals to zero, and letting the next guess for θ θ theta\theta 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 f ff plotted along with the line y = 0 y = 0 y=0y=0. We're trying to find θ θ theta\theta so that f ( θ ) = 0 f ( θ ) = 0 f(theta)=0f(\theta)=0; the value of θ θ theta\theta that achieves this is about 1.3 1.3 1.31.3. Suppose we initialized the algorithm with θ = 4.5 θ = 4.5 theta=4.5\theta=4.5. Newton's method then fits a straight line tangent to f f ff at θ = 4.5 θ = 4.5 theta=4.5\theta=4.5, and solves for the where that line evaluates to 0 0 00. (Middle figure.) This give us the next guess for θ θ theta\theta, which is about 2.8 2.8 2.82.8. The rightmost figure shows the result of running one more iteration, which the updates θ θ theta\theta to about 1.8 1.8 1.81.8. After a few more iterations, we rapidly approach θ = 1.3 θ = 1.3 theta=1.3\theta=1.3.
Newton's method gives a way of getting to f ( θ ) = 0 f ( θ ) = 0 f(theta)=0f(\theta)=0. What if we want to use it to maximize some function ℓ ℓ ℓ\ell? The maxima of ℓ ℓ ℓ\ell correspond to points where its first derivative ℓ ′ ( θ ) ℓ ′ ( θ ) ℓ^(')(theta)\ell^{\prime}(\theta) is zero. So, by letting f ( θ ) = ℓ ′ ( θ ) f ( θ ) = ℓ ′ ( θ ) f(theta)=ℓ^(')(theta)f(\theta)=\ell^{\prime}(\theta), we can use the same algorithm to maximize ℓ ℓ ℓ\ell, and we obtain update rule:
θ := θ − ℓ ′ ( θ ) ℓ ′ ′ ( θ ) . θ := θ − ℓ ′ ( θ ) ℓ ′ ′ ( θ ) . theta:=theta-(ℓ^(')(theta))/(ℓ^('')(theta)).\theta:=\theta-\frac{\ell^{\prime}(\theta)}{\ell^{\prime \prime}(\theta)} .
(Something to think about: How would this change if we wanted to use Newton's method to minimize rather than maximize a function?)
Lastly, in our logistic regression setting, θ θ theta\theta 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 ∇ θ ℓ ( θ ) . θ := θ − H − 1 ∇ θ ℓ ( θ ) . theta:=theta-H^(-1)grad_(theta)ℓ(theta).\theta:=\theta-H^{-1} \nabla_{\theta} \ell(\theta) .
Here, ∇ θ ℓ ( θ ) ∇ θ ℓ ( θ ) grad_(theta)ℓ(theta)\nabla_{\theta} \ell(\theta) is, as usual, the vector of partial derivatives of ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta) with respect to the θ i θ i theta_(i)\theta_{i}'s; and H H HH is an d d dd-by- d d dd matrix (actually, d + 1 − b y − d + 1 d + 1 − b y − d + 1 d+1-by-d+1d+1-by-\mathrm{d}+1, assuming that we include the intercept term) called the Hessian, whose entries are given by
H i j = ∂ 2 ℓ ( θ ) ∂ θ i ∂ θ j . H i j = ∂ 2 ℓ ( θ ) ∂ θ i ∂ θ j . H_(ij)=(del^(2)ℓ(theta))/(deltheta_(i)deltheta_(j)).H_{i j}=\frac{\partial^{2} \ell(\theta)}{\partial \theta_{i} \partial \theta_{j}} .
Newton's method typically enjoys faster convergence than (batch) gradient 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 d dd-by- d d dd Hessian; but so long as d d dd is not too large, it is usually much faster overall. When Newton's method is applied to maximize the logistic regression log likelihood function ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta), the resulting method is also called Fisher scoring.
Part III

Generalized Linear Models[5]

So far, we've seen a regression example, and a classification example. In the regression example, we had y | x ; θ ∼ N ( μ , σ 2 ) y | x ; θ ∼ N μ , σ 2 y|x;theta∼N(mu,sigma^(2))y | x ; \theta \sim \mathcal{N}\left(\mu, \sigma^{2}\right), and in the classification one, y | x ; θ ∼ Bernoulli ⁡ ( ϕ ) y | x ; θ ∼ Bernoulli ⁡ ( ϕ ) y|x;theta∼Bernoulli(phi)y | x ; \theta \sim \operatorname{Bernoulli}(\phi), for some appropriate definitions of μ μ mu\mu and ϕ ϕ phi\phi as functions of x x xx and θ θ theta\theta. 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). We will also show how other models in the GLM family can be derived and applied to other classification and regression problems.

8. 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
(3) p ( y ; η ) = b ( y ) exp ⁡ ( η T T ( y ) − a ( η ) ) (3) p ( y ; η ) = b ( y ) exp ⁡ ( η T T ( y ) − a ( η ) ) {:(3)p(y;eta)=b(y)exp(eta^(T)T(y)-a(eta)):}\begin{equation} p(y ; \eta)=b(y) \exp (\eta^{T} T(y)-a(\eta)) \tag{3} \end{equation}
Here, η η eta\eta is called the natural parameter (also called the canonical parameter) of the distribution; T ( y ) T ( y ) T(y)T(y) is the sufficient statistic (for the distributions we consider, it will often be the case that T ( y ) = y T ( y ) = y T(y)=yT(y)=y); and a ( η ) a ( η ) a(eta)a(\eta) is the log partition function. The quantity e − a ( η ) e − a ( η ) e^(-a(eta))e^{-a(\eta)} essentially plays the role of a normalization constant, that makes sure the distribution p ( y ; η ) p ( y ; η ) p(y;eta)p(y ; \eta) sums/integrates over y y yy to 1 1 11.
A fixed choice of T , a T , a T,aT, a and b b bb defines a family (or set) of distributions that is parameterized by η η eta\eta; as we vary η η eta\eta, we then get different distributions within this family.
We now show that the Bernoulli and the Gaussian distributions are examples of exponential family distributions. The Bernoulli distribution with mean ϕ ϕ phi\phi, written Bernoulli ( ϕ ) ( ϕ ) (phi)(\phi), specifies a distribution over y ∈ { 0 , 1 } y ∈ { 0 , 1 } y in{0,1}y \in\{0,1\}, so that p ( y = 1 ; ϕ ) = ϕ ; p ( y = 0 ; ϕ ) = 1 − ϕ p ( y = 1 ; ϕ ) = ϕ ; p ( y = 0 ; ϕ ) = 1 − ϕ p(y=1;phi)=phi;p(y=0;phi)=1-phip(y=1 ; \phi)=\phi ; p(y=0 ; \phi)=1-\phi. As we vary ϕ ϕ phi\phi, we obtain Bernoulli distributions with different means. We now show that this class of Bernoulli distributions, ones obtained by varying ϕ ϕ phi\phi, is in the exponential family; i.e., that there is a choice of T , a T , a T,aT, a and b b bb so that Equation (3) 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 − ϕ ) ) . p ( y ; ϕ ) = ϕ y ( 1 − ϕ ) 1 − y = exp ⁡ ( y log ⁡ ϕ + ( 1 − y ) log ⁡ ( 1 − ϕ ) ) = exp ⁡ log ⁡ ϕ 1 − ϕ y + log ⁡ ( 1 − ϕ ) . {:[p(y;phi)=phi^(y)(1-phi)^(1-y)],[=exp(y log phi+(1-y)log(1-phi))],[=exp((log((phi)/(1-phi)))y+log(1-phi)).]:}\begin{aligned} p(y ; \phi) &=\phi^{y}(1-\phi)^{1-y} \\ &=\exp (y \log \phi+(1-y) \log (1-\phi)) \\ &=\exp \left(\left(\log \left(\frac{\phi}{1-\phi}\right)\right) y+\log (1-\phi)\right) . \end{aligned}
Thus, the natural parameter is given by η = log ⁡ ( ϕ / ( 1 − ϕ ) ) η = log ⁡ ( ϕ / ( 1 − ϕ ) ) eta=log(phi//(1-phi))\eta=\log (\phi /(1-\phi)). Interestingly, if we invert this definition for η η eta\eta by solving for ϕ ϕ phi\phi in terms of η η eta\eta, we obtain ϕ = ϕ = phi=\phi= 1 / ( 1 + e − η ) 1 / 1 + e − η 1//(1+e^(-eta))1 /\left(1+e^{-\eta}\right). 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 T ( y ) = y a ( η ) = − log ⁡ ( 1 − ϕ ) = log ⁡ 1 + e η b ( y ) = 1 {:[T(y)=y],[a(eta)=-log(1-phi)],[=log(1+e^(eta))],[b(y)=1]:}\begin{aligned} T(y) &=y \\ a(\eta) &=-\log (1-\phi) \\ &=\log \left(1+e^{\eta}\right) \\ b(y) &=1 \end{aligned}
This shows that the Bernoulli distribution can be written in the form of Equation (3), using an appropriate choice of T , a T , a T,aT, a and b b bb.
Let's now move on to consider the Gaussian distribution. Recall that, when deriving linear regression, the value of σ 2 σ 2 sigma^(2)\sigma^{2} had no effect on our final choice of θ θ theta\theta and h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x). Thus, we can choose an arbitrary value for σ 2 σ 2 sigma^(2)\sigma^{2} without changing anything. To simplify the derivation below, let's set σ 2 = 1 σ 2 = 1 sigma^(2)=1\sigma^{2}=1.[6] We then have:
p ( y ; μ ) = 1 2 π exp ⁡ ( − 1 2 ( y − μ ) 2 ) = 1 2 π exp ⁡ ( − 1 2 y 2 ) ⋅ exp ⁡ ( μ y − 1 2 μ 2 ) p ( y ; μ ) = 1 2 π exp ⁡ − 1 2 ( y − μ ) 2 = 1 2 π exp ⁡ − 1 2 y 2 ⋅ exp ⁡ μ y − 1 2 μ 2 {:[p(y;mu)=(1)/(sqrt(2pi))exp(-(1)/(2)(y-mu)^(2))],[=(1)/(sqrt(2pi))exp(-(1)/(2)y^(2))*exp(mu y-(1)/(2)mu^(2))]:}\begin{aligned} p(y ; \mu) &=\frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{1}{2}(y-\mu)^{2}\right) \\ &=\frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{1}{2} y^{2}\right) \cdot \exp \left(\mu y-\frac{1}{2} \mu^{2}\right) \end{aligned}
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 ) . η = μ T ( y ) = y a ( η ) = μ 2 / 2 = η 2 / 2 b ( y ) = ( 1 / 2 π ) exp ⁡ − y 2 / 2 . {:[eta=mu],[T(y)=y],[a(eta)=mu^(2)//2],[=eta^(2)//2],[b(y)=(1//sqrt(2pi))exp(-y^(2)//2).]:}\begin{aligned} \eta &=\mu \\ T(y) &=y \\ a(\eta) &=\mu^{2} / 2 \\ &=\eta^{2} / 2 \\ b(y) &=(1 / \sqrt{2 \pi}) \exp \left(-y^{2} / 2\right) . \end{aligned}
There're many other distributions that are members of the exponential family: The multinomial (which we'll see later), the Poisson (for modelling count-data; also see the problem set); the gamma and the exponential (for modelling continuous, non-negative random variables, such as timeintervals); 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 y yy (given x x xx and θ θ theta\theta) comes from any of these distributions.

9. Constructing GLMs

Suppose you would like to build a model to estimate the number y y yy of customers arriving in your store (or number of page-views on your website) in any given hour, based on certain features x x xx such as store promotions, recent advertising, weather, day-of-week, etc. We know that the Poisson distribution 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 y yy as a function of x x xx. To derive a GLM for this problem, we will make the following three assumptions about the conditional distribution of y y yy given x x xx and about our model:
  1. y ∣ x ; θ ∼ ExponentialFamily ( η ) y ∣ x ; θ ∼ ExponentialFamily ( η ) y∣x;theta∼"ExponentialFamily"(eta)y \mid x ; \theta \sim\text{ExponentialFamily}(\eta). I.e., given x x xx and θ θ theta\theta, the distribution of y y yy follows some exponential family distribution, with parameter η η eta\eta.
  2. Given x x xx, our goal is to predict the expected value of T ( y ) T ( y ) T(y)T(y) given x x xx. In most of our examples, we will have T ( y ) = y T ( y ) = y T(y)=yT(y)=y, so this means we would like the prediction h ( x ) h ( x ) h(x)h(x) output by our learned hypothesis h h hh to satisfy h ( x ) = E [ y ∣ x ] h ( x ) = E [ y ∣ x ] h(x)=E[y∣x]h(x)=\mathrm{E}[y \mid x]. (Note that this assumption is satisfied in the choices for h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x) for both logistic regression and linear regression. For instance, in logistic regression, we had h θ ( x ) = p ( y = 1 ∣ x ; θ ) = 0 ⋅ p ( y = h θ ( x ) = p ( y = 1 ∣ x ; θ ) = 0 ⋅ p ( y = h_(theta)(x)=p(y=1∣x;theta)=0*p(y=h_{\theta}(x)=p(y=1 \mid x ; \theta)=0 \cdot p(y= 0 ∣ x ; θ ) + 1 ⋅ p ( y = 1 ∣ x ; θ ) = E [ y ∣ x ; θ ] . ) 0 ∣ x ; θ ) + 1 ⋅ p ( y = 1 ∣ x ; θ ) = E [ y ∣ x ; θ ] . ) 0∣x;theta)+1*p(y=1∣x;theta)=E[y∣x;theta].)0 \mid x ; \theta)+1 \cdot p(y=1 \mid x ; \theta)=\mathrm{E}[y \mid x ; \theta] .)
  3. The natural parameter η η eta\eta and the inputs x x xx are related linearly: η = θ T x η = θ T x eta=theta^(T)x\eta=\theta^{T} x. (Or, if η η eta\eta is vector-valued, then η i = θ i T x η i = θ i T x eta_(i)=theta_(i)^(T)x\eta_{i}=\theta_{i}^{T} 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 y yy; for example, we will shortly show that both logistic regression and ordinary least squares can both be derived as GLMs.

9.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 y yy (also called the response variable in GLM terminology) is continuous, and we model the conditional distribution of y y yy given x x xx as a Gaussian N ( μ , σ 2 ) N μ , σ 2 N(mu,sigma^(2))\mathcal{N}\left(\mu, \sigma^{2}\right). (Here, μ μ mu\mu may depend x x xx.) So, we let the ExponentialFamily ( η ) ExponentialFamily ( η ) "ExponentialFamily"(eta)\textit{ExponentialFamily}(\eta) distribution above be the Gaussian distribution. As we saw previously, in the formulation of the Gaussian as an exponential family distribution, we had μ = η μ = η mu=eta\mu=\eta. So, we have
h θ ( x ) = E [ y ∣ x ; θ ] = μ = η = θ T x . h θ ( x ) = E [ y ∣ x ; θ ] = μ = η = θ T x . {:[h_(theta)(x)=E[y∣x;theta]],[=mu],[=eta],[=theta^(T)x.]:}\begin{aligned} h_{\theta}(x) &=E[y \mid x ; \theta] \\ &=\mu \\ &=\eta \\ &=\theta^{T} x . \end{aligned}
The first equality follows from Assumption 2, above; the second equality follows from the fact that y ∣ x ; θ ∼ N ( μ , σ 2 ) y ∣ x ; θ ∼ N μ , σ 2 y∣x;theta∼N(mu,sigma^(2))y \mid x ; \theta \sim \mathcal{N}\left(\mu, \sigma^{2}\right), and so its expected value is given by μ μ mu\mu; the third equality follows from Assumption 1 (and our earlier derivation showing that μ = η μ = η mu=eta\mu=\eta in the formulation of the Gaussian as an exponential family distribution); and the last equality follows from Assumption 3.

9.2. Logistic Regression

We now consider logistic regression. Here we are interested in binary classification, so y ∈ { 0 , 1 } y ∈ { 0 , 1 } y in{0,1}y \in\{0,1\}. Given that y y yy is binary-valued, it therefore seems natural to choose the Bernoulli family of distributions to model the conditional distribution of y y yy given x x xx. In our formulation of the Bernoulli distribution as an exponential family distribution, we had ϕ = 1 / ( 1 + e − η ) ϕ = 1 / 1 + e − η phi=1//(1+e^(-eta))\phi=1 /\left(1+e^{-\eta}\right). Furthermore, note that if y | x ; θ ∼ Bernoulli ⁡ ( ϕ ) y | x ; θ ∼ Bernoulli ⁡ ( ϕ ) y|x;theta∼Bernoulli(phi)y | x ; \theta \sim \operatorname{Bernoulli}(\phi), then E [ y | x ; θ ] = ϕ E [ y | x ; θ ] = ϕ E[y|x;theta]=phi\mathrm{E}[y | x ; \theta]=\phi. 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 ) h θ ( x ) = E [ y ∣ x ; θ ] = ϕ = 1 / 1 + e − η = 1 / 1 + e − θ T x {:[h_(theta)(x)=E[y∣x;theta]],[=phi],[=1//(1+e^(-eta))],[=1//(1+e^(-theta^(T)x))]:}\begin{aligned} h_{\theta}(x) &=E[y \mid x ; \theta] \\ &=\phi \\ &=1 /\left(1+e^{-\eta}\right) \\ &=1 /\left(1+e^{-\theta^{T} x}\right) \end{aligned}
So, this gives us hypothesis functions of the form h θ ( x ) = 1 / ( 1 + e − θ T x ) h θ ( x ) = 1 / 1 + e − θ T x h_(theta)(x)=1//(1+e^(-theta^(T)x))h_{\theta}(x)=1 /\left(1+e^{-\theta^{T} x}\right). If you are previously wondering how we came up with the form of the logistic function 1 / ( 1 + e − z ) 1 / 1 + e − z 1//(1+e^(-z))1 /\left(1+e^{-z}\right), this gives one answer: Once we assume that y y yy conditioned on x x xx 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 g gg giving the distribution's mean as a function of the natural parameter ( g ( η ) = E [ T ( y ) ; η ] ) ( g ( η ) = E [ T ( y ) ; η ] ) (g(eta)=E[T(y);eta])(g(\eta)=\mathrm{E}[T(y) ; \eta]) is called the canonical response function. Its inverse, g − 1 g − 1 g^(-1)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.[7]

9.3. Softmax Regression

Let's look at one more example of a GLM. Consider a classification problem in which the response variable y y yy can take on any one of k k kk values, so y ∈ y ∈ y iny \in { 1 , 2 , … , k } { 1 , 2 , … , k } {1,2,dots,k}\{1,2, \ldots, k\}. For example, rather than classifying email into the two classes 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 k kk possible outcomes, one could use k k kk parameters ϕ 1 , … , ϕ k ϕ 1 , … , ϕ k phi_(1),dots,phi_(k)\phi_{1}, \ldots, \phi_{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 k − 1 k-1k-1 of the ϕ i ϕ i phi_(i)\phi_{i}'s uniquely determines the last one, as they must satisfy ∑ i = 1 k ϕ i = 1 ∑ i = 1 k   ϕ i = 1 sum_(i=1)^(k)phi_(i)=1\sum_{i=1}^{k} \phi_{i}=1). So, we will instead parameterize the multinomial with only k − 1 k − 1 k-1k-1 parameters, ϕ 1 , … , ϕ k − 1 ϕ 1 , … , ϕ k − 1 phi_(1),dots,phi_(k-1)\phi_{1}, \ldots, \phi_{k-1}, where ϕ i = p ( y = i ; ϕ ) ϕ i = p ( y = i ; ϕ ) phi_(i)=p(y=i;phi)\phi_{i}=p(y=i ; \phi), and p ( y = k ; ϕ ) = 1 − ∑ i = 1 k − 1 ϕ i p ( y = k ; ϕ ) = 1 − ∑ i = 1 k − 1   ϕ i p(y=k;phi)=1-sum_(i=1)^(k-1)phi_(i)p(y=k ; \phi)=1-\sum_{i=1}^{k-1} \phi_{i}. For notational convenience, we will also let ϕ k = 1 − ∑ i = 1 k − 1 ϕ i ϕ k = 1 − ∑ i = 1 k − 1   ϕ i phi_(k)=1-sum_(i=1)^(k-1)phi_(i)\phi_{k}=1-\sum_{i=1}^{k-1} \phi_{i}, but we should keep in mind that this is not a parameter, and that it is fully specified by ϕ 1 , … , ϕ k − 1 ϕ 1 , … , ϕ k − 1 phi_(1),dots,phi_(k-1)\phi_{1}, \ldots, \phi_{k-1}.
To express the multinomial as an exponential family distribution, we will define T ( y ) ∈ R k − 1 T ( y ) ∈ R k − 1 T(y)inR^(k-1)T(y) \in \mathbb{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 ] , 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 , T(1)=[[1],[0],[0],[vdots],[0]],T(2)=[[0],[1],[0],[vdots],[0]],T(3)=[[0],[0],[1],[vdots],[0]],cdots,T(k-1)=[[0],[0],[0],[vdots],[1]],T(k)=[[0],[0],[0],[vdots],[0]],T(1)=\left[\begin{array}{c}1 \\ 0 \\ 0 \\ \vdots \\ 0\end{array}\right], T(2)=\left[\begin{array}{c}0 \\ 1 \\ 0 \\ \vdots \\ 0\end{array}\right], T(3)=\left[\begin{array}{c}0 \\ 0 \\ 1 \\ \vdots \\ 0\end{array}\right], \cdots, T(k-1)=\left[\begin{array}{c}0 \\ 0 \\ 0 \\ \vdots \\ 1\end{array}\right], T(k)=\left[\begin{array}{c}0 \\ 0 \\ 0 \\ \vdots \\ 0\end{array}\right],
Unlike our previous examples, here we do not have T ( y ) = y T ( y ) = y T(y)=yT(y)=y; also, T ( y ) T ( y ) T(y)T(y) is now a k − 1 k − 1 k-1k-1 dimensional vector, rather than a real number. We will write ( T ( y ) ) i ( T ( y ) ) i (T(y))_(i)(T(y))_{i} to denote the i i ii-th element of the vector T ( y ) T ( y ) T(y)T(y).
We introduce one more very useful piece of notation. An indicator function 1 { ⋅ } 1 { ⋅ } 1{*}1\{\cdot\} takes on a value of 1 1 11 if its argument is true, and 0 0 00 otherwise ( 1 { ( 1 { (1{(1\{ True } = 1 , 1 { } = 1 , 1 { }=1,1{\}=1,1\{ False } = 0 ) } = 0 ) }=0)\}=0). For example, 1 { 2 = 3 } = 0 1 { 2 = 3 } = 0 1{2=3}=01\{2=3\}=0, and 1 { 3 = 1 { 3 = 1{3=1\{3= 5 − 2 } = 1 5 − 2 } = 1 5-2}=15-2\}=1. So, we can also write the relationship between T ( y ) T ( y ) T(y)T(y) and y y yy as ( T ( y ) ) i = 1 { y = i } ( T ( y ) ) i = 1 { y = i } (T(y))_(i)=1{y=i}(T(y))_{i}=1\{y=i\}. (Before you continue reading, please make sure you understand why this is true!) Further, we have that E [ ( T ( y ) ) i ] = P ( y = i ) = ϕ i E ( T ( y ) ) i = P ( y = i ) = ϕ i E[(T(y))_(i)]=P(y=i)=phi_(i)\mathrm{E}\left[(T(y))_{i}\right]=P(y=i)=\phi_{i}.
We are now ready to show that the multinomial is a member of the exponential family. We have:
p ( y ; ϕ ) = ϕ 1 1 { y = 1 } ϕ 2 1 { y = 2 } ⋯ ϕ k 1 { y = k } = ϕ 1 1 { y = 1 } ϕ 2 1 { y = 2 } ⋯ ϕ k 1 − ∑ i = 1 k − 1 1 { y = i } = ϕ 1 ( T ( y ) ) 1 ϕ 2 ( T ( y ) ) 2 ⋯ ϕ k 1 − ∑ i = 1 k − 1 ( T ( y ) ) i = exp ⁡ ( ( T ( y ) ) 1 log ⁡ ( ϕ 1 ) + ( T ( y ) ) 2 log ⁡ ( ϕ 2 ) + ⋯ + ( 1 − ∑ i = 1 k − 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 ( η ) ) p ( y ; ϕ ) = ϕ 1 1 { y = 1 } ϕ 2 1 { y = 2 } ⋯ ϕ k 1 { y = k } = ϕ 1 1 { y = 1 } ϕ 2 1 { y = 2 } ⋯ ϕ k 1 − ∑ i = 1 k − 1   1 { y = i } = ϕ 1 ( T ( y ) ) 1 ϕ 2 ( T ( y ) ) 2 ⋯ ϕ k 1 − ∑ i = 1 k − 1   ( T ( y ) ) i = exp ⁡ ( T ( y ) ) 1 log ⁡ ϕ 1 + ( T ( y ) ) 2 log ⁡ ϕ 2 + ⋯ + 1 − ∑ i = 1 k − 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 ( η ) ) {:[p(y;phi)=phi_(1)^(1{y=1})phi_(2)^(1{y=2})cdotsphi_(k)^(1{y=k})],[=phi_(1)^(1{y=1})phi_(2)^(1{y=2})cdotsphi_(k)^(1-sum_(i=1)^(k-1)1{y=i})],[=phi_(1)^((T(y))_(1))phi_(2)^((T(y))_(2))cdotsphi_(k)^(1-sum_(i=1)^(k-1)(T(y))_(i))],[=exp((T(y))_(1)log(phi_(1))+(T(y))_(2)log(phi_(2))+:}],[qquad cdots+(1-sum_(i=1)^(k-1)(T(y))_(i))log(phi_(k)))],[=exp((T(y))_(1)log(phi_(1)//phi_(k))+(T(y))_(2)log(phi_(2)//phi_(k))+],[qquad cdots+(T(y))_(k-1)log(phi_(k-1)//phi_(k))+log(phi_(k)))],[=b(y)exp(eta^(T)T(y)-a(eta))]:}\begin{aligned} p(y ; \phi)=& \phi_{1}^{1\{y=1\}} \phi_{2}^{1\{y=2\}} \cdots \phi_{k}^{1\{y=k\}} \\ =& \phi_{1}^{1\{y=1\}} \phi_{2}^{1\{y=2\}} \cdots \phi_{k}^{1-\sum_{i=1}^{k-1} 1\{y=i\}} \\ =& \phi_{1}^{(T(y))_{1}} \phi_{2}^{(T(y))_{2}} \cdots \phi_{k}^{1-\sum_{i=1}^{k-1}(T(y))_{i}} \\ =& \exp \left((T(y))_{1} \log \left(\phi_{1}\right)+(T(y))_{2} \log \left(\phi_{2}\right)+\right.\\ &\qquad \cdots+\left(1-\sum_{i=1}^{k-1}(T(y))_{i}\right) \log (\phi_{k})) \\ =& \exp ((T(y))_{1} \log (\phi_{1} / \phi_{k})+(T(y))_{2} \log (\phi_{2} / \phi_{k})+\\ &\qquad \cdots+(T(y))_{k-1} \log (\phi_{k-1} / \phi_{k})+\log (\phi_{k})) \\ =& b(y) \exp (\eta^{T} T(y)-a(\eta)) \end{aligned}
where
η = [ log ⁡ ( ϕ 1 / ϕ k ) log ⁡ ( ϕ 2 / ϕ k ) ⋮ log ⁡ ( ϕ k − 1 / ϕ k ) ] , a ( η ) = − log ⁡ ( ϕ k ) b ( y ) = 1 . η = log ⁡ ϕ 1 / ϕ k log ⁡ ϕ 2 / ϕ k ⋮ log ⁡ ϕ k − 1 / ϕ k , a ( η ) = − log ⁡ ϕ k b ( y ) = 1 . {:[eta=[[log(phi_(1)//phi_(k))],[log(phi_(2)//phi_(k))],[vdots],[log(phi_(k-1)//phi_(k))]]","],[a(eta)=-log(phi_(k))],[b(y)=1.]:}\begin{aligned} \eta &=\left[\begin{array}{c} \log \left(\phi_{1} / \phi_{k}\right) \\ \log \left(\phi_{2} / \phi_{k}\right) \\ \vdots \\ \log \left(\phi_{k-1} / \phi_{k}\right) \end{array}\right], \\ a(\eta) &=-\log \left(\phi_{k}\right) \\ b(y) &=1 . \end{aligned}
This completes our formulation of the multinomial as an exponential family distribution.
The link function is given (for i = 1 , … , k i = 1 , … , k i=1,dots,ki=1, \ldots, k) by
η i = log ⁡ ϕ i ϕ k . η i = log ⁡ ϕ i ϕ k . eta_(i)=log (phi_(i))/(phi_(k)).\eta_{i}=\log \frac{\phi_{i}}{\phi_{k}} .
For convenience, we have also defined η k = log ⁡ ( ϕ k / ϕ k ) = 0 η k = log ⁡ ϕ k / ϕ k = 0 eta_(k)=log(phi_(k)//phi_(k))=0\eta_{k}=\log \left(\phi_{k} / \phi_{k}\right)=0. To invert the link function and derive the response function, we therefore have that
(4) e η i = ϕ i ϕ k ϕ k e η i = ϕ i ϕ k ∑ i = 1 k e η i = ∑ i = 1 k ϕ i = 1 (4) e η i = ϕ i ϕ k ϕ k e η i = ϕ i ϕ k ∑ i = 1 k   e η i = ∑ i = 1 k   ϕ i = 1 {:(4){:[e^(eta_(i))=(phi_(i))/(phi_(k))],[phi_(k)e^(eta_(i))=phi_(i)],[phi_(k)sum_(i=1)^(k)e^(eta_(i))=sum_(i=1)^(k)phi_(i)=1]:}:}\begin{equation} \begin{split} e^{\eta_{i}} &=\frac{\phi_{i}}{\phi_{k}} \\ \phi_{k} e^{\eta_{i}} &=\phi_{i} \\ \phi_{k} \sum_{i=1}^{k} e^{\eta_{i}} &=\sum_{i=1}^{k} \phi_{i}=1 \end{split} \tag{4} \end{equation}
This implies that ϕ k = 1 / ∑ i = 1 k e η i ϕ k = 1 / ∑ i = 1 k   e η i phi_(k)=1//sum_(i=1)^(k)e^(eta_(i))\phi_{k}=1 / \sum_{i=1}^{k} e^{\eta_{i}}, which can be substituted back into Equation (4) to give the response function
ϕ i = e η i ∑ j = 1 k e η j ϕ i = e η i ∑ j = 1 k   e η j phi_(i)=(e^(eta_(i)))/(sum_(j=1)^(k)e^(eta_(j)))\phi_{i}=\frac{e^{\eta_{i}}}{\sum_{j=1}^{k} e^{\eta_{j}}}
This function mapping from the η η eta\eta's to the ϕ ϕ phi\phi's is called the softmax function.
To complete our model, we use Assumption 3, given earlier, that the η i η i eta_(i)\eta_{i}'s are linearly related to the x x xx's. So, have η i = θ i T x η i = θ i T x eta_(i)=theta_(i)^(T)x\eta_{i}=\theta_{i}^{T} x (for i = 1 , … , k − 1 i = 1 , … , k − 1 i=1,dots,k-1i=1, \ldots, k-1 ), where θ 1 , … , θ k − 1 ∈ R d + 1 θ 1 , … , θ k − 1 ∈ R d + 1 theta_(1),dots,theta_(k-1)inR^(d+1)\theta_{1}, \ldots, \theta_{k-1} \in \mathbb{R}^{d+1} are the parameters of our model. For notational convenience, we can also define θ k = 0 θ k = 0 theta_(k)=0\theta_{k}=0, so that η k = θ k T x = 0 η k = θ k T x = 0 eta_(k)=theta_(k)^(T)x=0\eta_{k}=\theta_{k}^{T} x=0, as given previously. Hence, our model assumes that the conditional distribution of y y yy given x x xx is given by
p ( y = i | x ; θ ) = ϕ i = e η i ∑ j = 1 k e η j = (5) e θ i T x ∑ j = 1 k e θ j T x p ( y = i | x ; θ )      = ϕ i      = e η i ∑ j = 1 k   e η j      = (5) e θ i T x ∑ j = 1 k   e θ j T x {:[p(y=i|x;theta)=phi_(i)],[=(e^(eta_(i)))/(sum_(j=1)^(k)e^(eta_(j)))],[={:(5)(e^(theta_(i)^(T)x))/(sum_(j=1)^(k)e^(theta_(j)^(T)x)):}]:}\begin{split} p(y=i | x ; \theta) &=\phi_{i} \\ &=\frac{e^{\eta_{i}}}{\sum_{j=1}^{k} e^{\eta_{j}}} \\ &=\begin{equation}\frac{e^{\theta_{i}^{T} x}}{\sum_{j=1}^{k} e^{\theta_{j}^{T} x}}\tag{5}\end{equation} \end{split}
This model, which applies to classification problems where y ∈ { 1 , … , k } y ∈ { 1 , … , k } y in{1,dots,k}y \in\{1, \ldots, 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 ⁡ ( θ 1 T x ) ∑ j = 1 k exp ⁡ ( θ j T x ) exp ⁡ ( θ 2 T x ) ∑ j = 1 k exp ⁡ ( θ j T x ) ⋮ exp ⁡ ( θ k − 1 T x ) ∑ j = 1 k exp ⁡ ( θ j T x ) ] . h θ ( x ) = E [ T ( y ) | x ; θ ] = E 1 { y = 1 } 1 { y = 2 } ⋮ 1 { y = k − 1 } x ; θ = ϕ 1 ϕ 2 ⋮ ϕ k − 1 = exp ⁡ θ 1 T x ∑ j = 1 k   exp ⁡ θ j T x exp ⁡ θ 2 T x ∑ j = 1 k   exp ⁡ θ j T x ⋮ exp ⁡ θ k − 1 T x ∑ j = 1 k   exp ⁡ θ j T x . {:[h_(theta)(x)=E[T(y)|x;theta]],[=E[[1{y=1}],[1{y=2}],[vdots],[1{y=k-1}]|x;theta]],[=[[phi_(1)],[phi_(2)],[vdots],[phi_(k-1)]]],[=[[(exp(theta_(1)^(T)x))/(sum_(j=1)^(k)exp(theta_(j)^(T)x))],[(exp(theta_(2)^(T)x))/(sum_(j=1)^(k)exp(theta_(j)^(T)x))],[vdots],[(exp(theta_(k-1)^(T)x))/(sum_(j=1)^(k)exp(theta_(j)^(T)x))]].]:}\begin{aligned} h_{\theta}(x)&= \mathrm{E}[T(y) |x ; \theta] \\ &= \mathrm{E}\left[\left.\begin{array}{c} 1\{y=1\} \\ 1\{y=2\} \\ \vdots \\ 1\{y=k-1\} \end{array} \right| x ; \theta\right] \\ &= \left[\begin{array}{c} \phi_{1} \\ \phi_{2} \\ \vdots \\ \phi_{k-1} \end{array}\right] \\ &= \left[\begin{array}{c} \frac{\exp \left(\theta_{1}^{T} x\right)}{\sum_{j=1}^{k} \exp \left(\theta_{j}^{T} x\right)} \\ \frac{\exp \left(\theta_{2}^{T} x\right)}{\sum_{j=1}^{k} \exp \left(\theta_{j}^{T} x\right)} \\ \vdots \\ \frac{\exp \left(\theta_{k-1}^{T} x\right)}{\sum_{j=1}^{k} \exp \left(\theta_{j}^{T} x\right)} \end{array}\right]. \end{aligned}
In other words, our hypothesis will output the estimated probability that p ( y = i ∣ x ; θ ) p ( y = i ∣ x ; θ ) p(y=i∣x;theta)p(y=i \mid x ; \theta), for every value of i = 1 , … , k i = 1 , … , k i=1,dots,ki=1, \ldots, k. (Even though h θ ( x ) h θ ( x ) h_(theta)(x)h_{\theta}(x) as defined above is only k − 1 k − 1 k-1k-1 dimensional, clearly p ( y = k ∣ x ; θ ) p ( y = k ∣ x ; θ ) p(y=k∣x;theta)p(y=k \mid x ; \theta) can be obtained as 1 − ∑ i = 1 k − 1 ϕ i . ) 1 − ∑ i = 1 k − 1   ϕ i . {:1-sum_(i=1)^(k-1)phi_(i).)\left.1-\sum_{i=1}^{k-1} \phi_{i} .\right)
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 n nn examples { ( x ( i ) , y ( i ) ) ; i = 1 , … , n } x ( i ) , y ( i ) ; i = 1 , … , n {(x^((i)),y^((i)));i=1,dots,n}\left\{\left(x^{(i)}, y^{(i)}\right) ; i=1, \ldots, n\right\} and would like to learn the parameters θ i θ i theta_(i)\theta_{i} of this model, we would begin by writing down the log-likelihood
ℓ ( θ ) = ∑ i = 1 n log ⁡ p ( y ( i ) ∣ x ( i ) ; θ ) = ∑ i = 1 n log ⁡ ∏ l = 1 k ( e l θ l T x ( i ) ∑ j = 1 k e θ j T x ( i ) ) 1 { y ( i ) = l } ℓ ( θ ) = ∑ i = 1 n   log ⁡ p y ( i ) ∣ x ( i ) ; θ = ∑ i = 1 n   log ⁡ ∏ l = 1 k   e l θ l T x ( i ) ∑ j = 1 k   e θ j T x ( i ) 1 y ( i ) = l {:[ℓ(theta)=sum_(i=1)^(n)log p(y^((i))∣x^((i));theta)],[=sum_(i=1)^(n)log prod_(l=1)^(k)((e_(l)^(theta_(l)^(T)x^((i))))/(sum_(j=1)^(k)e^(theta_(j)^(T)x^((i)))))^(1{y^((i))=l})]:}\begin{aligned} \ell(\theta) &=\sum_{i=1}^{n} \log p\left(y^{(i)} \mid x^{(i)} ; \theta\right) \\ &=\sum_{i=1}^{n} \log \prod_{l=1}^{k}\left(\frac{e_{l}^{\theta_{l}^{T} x^{(i)}}}{\sum_{j=1}^{k} e^{\theta_{j}^{T} x^{(i)}}}\right)^{1\left\{y^{(i)}=l\right\}} \end{aligned}
To obtain the second line above, we used the definition for p ( y ∣ x ; θ ) p ( y ∣ x ; θ ) p(y∣x;theta)p(y \mid x ; \theta) given in Equation (5). We can now obtain the maximum likelihood estimate of the parameters by maximizing ℓ ( θ ) ℓ ( θ ) ℓ(theta)\ell(\theta) in terms of θ θ theta\theta, using a method such as gradient ascent or Newton's method.
This was the first lecture from Andrew Ng's CS229 course on Machine Learning. You can read the notes from the next lecture on Generative Learning Algorithms here.
  1. We use the notation " a := b a := b a:=ba:=b" to denote an operation (in a computer program) in which we set the value of a variable a a aa to be equal to the value of b b bb. In other words, this operation overwrites a a aa with the value of b b bb. In contrast, we will write " a = b a = b a=ba=b" when we are asserting a statement of fact, that the value of a a aa is equal to the value of b b bb. ↩︎
  2. By slowly letting the learning rate α α alpha\alpha decrease to zero as the algorithm runs, it is also possible to ensure that the parameters will converge to the global minimum rather than merely oscillate around the minimum. ↩︎
  3. Note that in the above step, we are implicitly assuming that X T X X T X X^(T)XX^{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 X T X X^(T)XX^{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. ↩︎
  4. If x x xx is vector-valued, this is generalized to be w ( i ) = exp ⁡ ( − ( x ( i ) − x ) T ( x ( i ) − x ) / ( 2 τ 2 ) ) w ( i ) = exp ⁡ ( − ( x ( i ) − x ) T ( x ( i ) − x ) / ( 2 τ 2 ) ) w^((i))=exp(-(x^((i))-x)^(T)(x^((i))-x)//(2tau^(2)))w^{(i)}=\exp (-(x^{(i)}-x)^{T}(x^{(i)}-x) /(2 \tau^{2})), or w ( i ) = exp ⁡ ( − ( x ( i ) − x ) T Σ − 1 ( x ( i ) − x ) / ( 2 τ 2 ) ) w ( i ) = exp ⁡ ( − ( x ( i ) − x ) T Σ − 1 ( x ( i ) − x ) / ( 2 τ 2 ) ) w^((i))=exp(-(x^((i))-x)^(T)Sigma^(-1)(x^((i))-x)//(2tau^(2)))w^{(i)}=\exp (-(x^{(i)}-x)^{T} \Sigma^{-1}(x^{(i)}-x) /(2 \tau^{2})), for an appropriate choice of τ τ tau\tau or Σ Σ Sigma\Sigma. ↩︎
  5. 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.). ↩︎
  6. If we leave σ 2 σ 2 sigma^(2)\sigma^{2} as a variable, the Gaussian distribution can also be shown to be in the exponential family, where η ∈ R 2 η ∈ R 2 eta inR^(2)\eta \in \mathbb{R}^{2} is now a 2 2 22-dimension vector that depends on both μ μ mu\mu and σ σ sigma\sigma. For the purposes of GLMs, however, the σ 2 σ 2 sigma^(2)\sigma^{2} parameter can also be treated by considering a more general definition of the exponential family: p ( y ; η , τ ) = b ( a , τ ) exp ⁡ ( ( η T T ( y ) − p ( y ; η , τ ) = b ( a , τ ) exp ⁡ η T T ( y ) − p(y;eta,tau)=b(a,tau)exp((eta^(T)T(y)-:}p(y ; \eta, \tau)=b(a, \tau) \exp \left(\left(\eta^{T} T(y)-\right.\right. a ( η ) ) / c ( τ ) ) a ( η ) ) / c ( τ ) ) a(eta))//c(tau))a(\eta)) / c(\tau)). Here, τ τ tau\tau is called the dispersion parameter, and for the Gaussian, c ( τ ) = σ 2 c ( τ ) = σ 2 c(tau)=sigma^(2)c(\tau)=\sigma^{2}; but given our simplification above, we won't need the more general definition for the examples we will consider here. ↩︎
  7. Many texts use g g gg to denote the link function, and g − 1 g − 1 g^(-1)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. ↩︎

Recommended for you

Emil Junker
Manipulative Attacks in Group Identification
Manipulative Attacks in Group Identification
This review provides an introduction to the group identification problem and gives an overview of the feasibility and computational complexity of manipulative attacks in group identification.
2 points
0 issues
353 views
1 points
0 issues
112 views
Contents