Machine Learning 1. Lecture 6. Linear Regression Models

What this lecture is about

We continue studying linear models, and today it is about regression and the method of least squares, which, generally speaking, we have been talking about since the first lecture. The plan is this. We start with the classic, textbook problem of multiple linear regression and see how it is solved with the singular value decomposition. We talk about regularization again. Then about principal component analysis: it too is an application of least squares, but to a quite different problem — not building a predictive model but transforming the feature space. And at the end we will bring everything together, and I will show a very unexpected insight about least squares that is rarely written or talked about.

Multiple linear regression

Let us start with the familiar notation.

Today the answers are real numbers, $Y=\mathbb{R}$. In least squares the difference $a(x_i,w)-y_i$ is the model's error on the object; for mathematical convenience we square it and sum over all the objects of the training sample. We could also assign weights to the objects — by default all objects are equally important — and we have already seen one trick with weights: the Nadaraya–Watson formula was also least squares, and the game there was precisely with the weights. Today we leave them alone; the weights are ones. And a terminological detail: the value of the functional at the minimum point is called the residual sum of squares (RSS); if you meet it in the literature, do not be alarmed.

Today the model is linear. Objects are described by $n$ numerical features, we simply sum them with weights, and the weights are the parameters of the model. Recall the applied problem from the first lecture: what makes up the price of a flat? Floor area, a balcony, a rubbish chute, the number of storeys, the number of rooms and many other collected factors. «Makes up» means, apparently, a linear combination of features. Probably not bad, though the model will later have to be generalized — we will consider the nonlinear generalization in the next lecture. But today, the classics.

The «objects–features» matrix: the rows are objects, $\ell$ of them, the columns are the features $f_1,\dots ,f_n$; $y$ is a column vector of length $\ell$, and the parameter vector $w$ is also a column. Moving to matrix notation is very convenient. $Fw$ is a column vector whose elements are the values of the model $a(x_i,w)$ on all the objects of the training sample, and the sum of squared differences is the squared norm of the vector $Fw-y$. We get a compact notation, and the way to minimize over $w$ suggests itself: differentiate and set the derivative to zero.

The normal system of equations

Differentiating with respect to a vector is not a very complicated thing; I will not go into its subtleties; you should have been told them in other courses.

whence follows the normal system of the least squares problem:

where $F^{\mathsf{T}}F$ is a symmetric $n\times n$ matrix. This is a system of $n$ linear equations in $n$ unknowns, and if the matrix $F^{\mathsf{T}}F$ is non-singular, the solution is written directly as

For the inverse matrix to exist, restrictions are needed, and we see them right from the course of the solution. There must be more objects than features, $\ell ⩾n$: otherwise the rank of the matrix will be less than $n$, and there will be no inverse. But even this does not guarantee the existence of the inverse: the features must be linearly independent, that is, $F$ must be a full-rank matrix, of rank equal to the number of its columns $n$.

Why «pseudoinverse»? Imagine that $F$ is square and non-singular, $n\times n$: as many objects as features. In our problem such a restriction has no meaningful sense — an extremely rare, unlikely situation that nobody needs; usually there are many more objects than features. But if $F$ were invertible after all, the brackets would open, $(AB{)}^{-1}=B^{-1}{A}^{-1}$, everything would cancel, and $F^{+}$ would turn out to be simply the inverse matrix. Then an exact solution would exist, $Fw=y$, with zero residual sum of squares, and we would find it by the same formula. So the pseudoinverse is a generalization of the inverse to the case where we invert not a square but, in general, a rectangular $\ell \times n$ matrix of full rank; if $F$ is square, it coincides with the ordinary inverse. And if we substitute $w^{\ast }$ into the functional, we get the matrix $P_F$ made of four matrices, the projection one: the vector $y$ is transformed somehow in the space ${\mathbb{R}}^{\ell }$ — of dimension equal to the number of objects — into its approximation $P_{F}y$. Let us work out this geometry.

The geometric meaning of least squares

The column $f_j\in {\mathbb{R}}^{\ell }$ of the matrix $F$ is the values of the $j$th feature on all the objects of the training sample. The linear span of the columns of the matrix $F=(f_1,\dots ,f_n)$:

is in effect the set of all possible linear models, a subspace in ${\mathbb{R}}^{\ell }$. $P_{F}y$ is the projection of the vector $y$ onto the subspace $\mathcal{L}(F)$, and $(I_{\ell }-P_F)y$ is the projection of $y$ onto its orthogonal complement.

Solving the least squares problem, we drop a perpendicular from the point $y$ in the space ${\mathbb{R}}^{\ell }$ onto the linear span of the columns of the matrix $F$. Least squares says: find in this subspace the point closest to $y$ — that is the projection of $y$ onto the linear span, and the matrix that transforms $y$ into this projection drops the perpendicular; that is why it is called the projection matrix. $P_{F}y$ is the vector of the model's answers on all the objects of the training sample, approximating the given $y$. The orthogonal projection $y-P_{F}y$ onto the orthogonal complement is the residual, the vector of the model's errors on all the objects. Clearly the original point would hardly have landed exactly in this subspace: the solution is always some approximation. This is important to understand.

The solution through the singular value decomposition

How is the pseudoinverse found? In different ways: there are a huge number of numerical methods implemented long ago, thousands of times, and telling how they work is outside our course. But I will drag you through one thing, because it is extremely useful, instructive and will come in handy in many other places.

The matrices $F{F}^{\mathsf{T}}$, of size $\ell \times \ell$, and $F^{\mathsf{T}}F$, of size $n\times n$, are very closely related: they have the same eigenvalues. More precisely, we agreed that $\ell >n$, both matrices have rank $n$, so the larger matrix $F{F}^{\mathsf{T}}$ has its remaining $\ell -n$ eigenvalues equal to zero, and the nonzero ones of the two matrices coincide. On the diagonal of $D$ are the roots of these common eigenvalues — the singular values of the matrix $F$. You have surely covered the singular value decomposition in courses on linear algebra and numerical methods, so this is a repetition.

Let us substitute this representation into the pseudoinverse. Remember that $U$ and $V$ are orthogonal: $V^{\mathsf{T}}V$ is the identity matrix and goes away at once; when the product is inverted, the $U$ cancel, $D^{-2}$ times $D$ remains, and they cancel too. Moving to writing in terms of the columns of the matrices $U$ and $V$, we get:

Identities: $(AB{)}^{-1}=B^{-1}{A}^{-1}$, $(AB{)}^{\mathsf{T}}=B^{\mathsf{T}}{A}^{\mathsf{T}}$, $\Vert w{\Vert }^2=w^{\mathsf{T}}w$.

$u_j{v}_j^{\mathsf{T}}$ is a column times a row, a rank-one matrix of the same size as $F^{+}$. The pseudoinverse has been represented as a sum of rank-one matrices divided by the singular values. Once the singular value decomposition has been computed, finding the pseudoinverse is no trouble at all. The vector of optimal parameters $w^{\ast }$ has been expanded in the basis of eigenvectors $u_j$ of the matrix $F^{\mathsf{T}}F$, and again ${\lambda }_j$ is in the denominator. The model's answers on the training sample, that is, the orthogonal projection onto the linear span of the columns of $F$, have been expanded in the eigenvectors of the matrix $F{F}^{\mathsf{T}}$ — and note: this is the only one of the four expressions without ${\lambda }_j$ in the denominator. The norm $\Vert w^{\ast }{\Vert }^2$ is again a sum, with ${\lambda }_j$ in the denominator.

Multicollinearity and the condition number

Let us reason about ${\lambda }_j$ in the denominator. What is bad about it? Recall the notion of an ill-conditioned matrix — a matrix that has eigenvalues close to zero. If some ${\lambda }_j$ tends to zero, a large unstable term appears in all the sums but one. You will see the effects of instability when you compute the pseudoinverse, when you compute the parameter vector $w^{\ast }$, and when you compute its norm: it will tend to plus infinity, which is a sign of an unstable solution. But you will see nothing of the kind when you look at how the model approximates the sample: nothing tends to any infinity there, and it will seem to you that the solution is great. Nevertheless it is unstable.

A singular value tending to zero says that there is an almost linear dependence between the columns of the matrix, that is, between the features. Were the dependence exact, the eigenvalue would be exactly zero, and we could not invert the matrix at all. Still more insidious is the situation where we can invert it, but the inversion is unstable because of ${\lambda }_j$ tending to zero in the denominator.

The condition number can be defined in different ways, proving that the definitions are equivalent; in particular, it is the ratio of the maximum eigenvalue to the minimum one. In our case $S=F^{\mathsf{T}}F$ — it is the one we invert — and $u=F^{\mathsf{T}}y$. Errors in measuring the features $f_j(x_i)$ and the answers $y_i$ are amplified by a factor of $\mu (F^{\mathsf{T}}F)$ when computing the least squares solution vector $w^{\ast }=(F^{\mathsf{T}}F{)}^{-1}{F}^{\mathsf{T}}y$. Let the relative measurement error be 1% and the condition number 1000, which may easily happen. Then the relative error of the model's parameter vector may reach 1000% — ten times the quantity itself. Having computed the singular value decomposition, we will see this with our own eyes: sort the eigenvalues in decreasing order, compare the minimum with the maximum, and if the condition number is large, something has to be done.

The problem of multicollinearity is closely tied to overfitting. If the matrix $F^{\mathsf{T}}F$ is ill-conditioned, the solution is unstable and poorly interpretable: the norm of the vector $w^{\ast }$ starts to tend to infinity, and we see coefficients large in absolute value and of different signs — a sure sign that the features are almost linearly dependent. On training it seems to us that everything is fine, but as soon as we use the model on an independent test sample $X^k$, the instability shows itself at once.

The first strategy is treatment that is, as doctors say, symptomatic: if the norm of the vector $w$ runs off to infinity, we do not let it, and from all solutions we choose the solution of minimum norm, the most stable one. The second: since the features are multicollinear and there are features depending on others, let us find independent ones and throw out the rest. The third: move to a smaller number of new features; we will consider it in the last part of the lecture. We solved the least squares problem in the first ten minutes, but considering the properties of the solution we see problems that have to be solved, including by changing the problem statement itself. Today we will look at all three strategies.

Ridge regression

The first strategy: penalize the model for increasing the $L_2$ norm of the weight vector. Among the $w$ that give the main functional values close to the minimum there may be too many solutions, and we will choose among them the solution of minimum norm: we add the two criteria with a regularization parameter.

The regularization parameter is sometimes conveniently written as $\tau$, a factor in front of the regularizer, and sometimes as $\frac{1}{\sigma }$. We will discuss probabilistic interpretations in the next lecture; for now I will only hint that there will be a Gaussian prior on $w$ here, and $\sigma$ is in effect its variance. Pulling the variance towards zero is the same as greatly increasing the regularization. Regularization is of course a double-edged sword: too small a $\tau$ does not solve the problem, too large a one drives the solution to zero. Differentiating, we get the addition $2\tau w$, and it results in $\tau I_n$ being added to the matrix $F^{\mathsf{T}}F$ — a «ridge» along the diagonal. Hence the name: ridge regression.

How will $\tau$ show up in the formulas through the singular value decomposition? Let us repeat everything we did before, with the addition:

The regularization coefficient has been added to ${\lambda }_j$ in the denominator and has thereby moved the denominator away from zero: if ${\lambda }_j$ tended to zero, we added $\tau$ to it and solved the problem. And the structure of the solution has stayed exactly the same. Let us check ourselves: with $\tau =0$ we get the same solution as before. With $\tau >0$, $F{w}_{\tau }^{\ast }\ne F{w}^{\ast }$, but the solution becomes much more stable. It becomes biased. The properties of bias are discussed in courses of mathematical statistics, where multiple linear regression is covered: the solution of the least squares problem is unbiased but may turn out unstable, and by adding regularization we traded unbiasedness for stability. As a result of this trade we get a more profitable solution, because stability does matter.

Choosing the regularization parameter

What is the advantage of the singular value decomposition? It can be computed once for the matrix $F$, giving the matrices $U$, $D$ and $V$ — that is a lot of computation, 99% of all the work — and then we can compute $w^{\ast }$, $F{w}^{\ast }$ and anything we want very quickly by the formulas, varying $\tau$ at that: the singular value decomposition does not depend on $\tau$. Too small a $\tau$ — weak regularization, instability; too large — a degenerate, zero solution. A compromise is needed, and we know how to tune such hyperparameters: let us do hold-out.

At the start the whole sample has to be split into training and control: the training one to build the model, and the control one, not very large, to optimize the regularization coefficient. We compute the criterion on the objects of the control sample, substituting the regularized solution, and get a remarkable formula that gives a very fast way to compute the dependence of the hold-out criterion on $\tau$. The singular value decomposition is computed once; $F^{\prime }U$ is multiplied once, $V^{\mathsf{T}}y$ also once. Now, having this $k\times n$ matrix and this vector, one can easily and quickly compute the criterion for any $\tau$: take $\tau$ on a grid, plot $Q(\tau )$, see the characteristic minimum, choose the optimal $\tau$. That is the charm of the singular value decomposition: having computed it once, we learn a great deal about the properties of the solution and get the ability to optimize the regularization coefficient.

Regularization reduces the «effective dimension»

Terminology: such regularization is also called shrinkage or weight decay:

People also say that a reduction of the effective dimension takes place in it. But we did not seem to reduce any dimension: we used $n$ features and keep using them. Why do people say this? There is an amusing fact: when a projection matrix projects a vector onto a subspace of smaller dimension, the trace of the projection matrix is the dimension of this subspace. The proof is trivial: under the trace sign matrices can be permuted, $\mathrm{tr}(AB)=\mathrm{tr}(BA)$, after which $F^{\mathsf{T}}F$ cancels with its inverse, and the trace of the identity matrix remains:

Let us do the same with regularization. We get a similar construction with the addition $\tau I_n$; it is not a projection matrix in the strict sense, it does not do an orthogonal projection, but it still looks for a vector in the space of linear combinations of the columns of $F$:

With $\tau =0$ we would again get $n$, while $\tau >0$ leads to a decrease in the trace. The trace here is not a dimension — it is not even an integer — but it is a kind of analogue of one, which is why people say the effective dimension decreases. Regularization squeezes the solution into a neighbourhood of zero: the norm of the weight vector shrinks, and a kind of analogue of the dimension decreases too. An amusing fact, needed just to understand the geometry of what is happening.

LASSO and non-smooth regularizers

So far we added the sum of squares $w_j^2$ to the functional. Can we add a sum of absolute values? We can. We have already discussed such a regularizer when we talked about its probabilistic interpretation, and today we have a regression problem.

This is LASSO, Least Absolute Shrinkage and Selection Operator. The name, of course, is fitted to the acronym, but the metaphor is rich and memorable: it is as if we throw a noose over the coefficient vector, and it is forced to get rid of excess degrees of freedom. Why can minimization with such a regularizer lead to some $w_j$ becoming zero? And a zero $w_j$ in a linear model means that the feature is always multiplied by zero and is no longer used — this is feature selection, that very second strategy.

For the unconstrained minimization problem with an $L_1$ regularizer there are methods of its own, but there is also a technique that reduces it to a smooth one. The first step: this statement is equivalent to a statement with an inequality constraint, where the regularization parameter is replaced by an upper bound $\varkappa$:

There is a monotonic dependence between $\mu$ and $\varkappa$; which exactly we do not discuss now, but the main thing is that for a given $\varkappa$ one can compute a $\mu$ for which the problems are equivalent. I will not prove it; the idea is this: write out the Karush–Kuhn–Tucker conditions for the problem with an inequality constraint (recall the last lecture) and make sure the resulting system is equivalent to the system of conditions for the original problem.

The second step: make the change of variables $w_j=w_j^{+}-w_j^{-}$, where both terms are non-negative and one of them is always zero. If $w_j>0$, then $w_j^{+}=w_j$ and $w_j^{-}=0$; if $w_j<0$, the other way round. Then the absolute value is expressed through the sum, $|w_j|=w_j^{+}+w_j^{-}$, and we get a smooth problem:

It can be solved with the same Karush–Kuhn–Tucker conditions, and complementary slackness conditions appear: the inequality constraints hold as equalities, they rest against their bounds. If we increase the regularization coefficient $\mu$, more and more constraints $w_j^{\pm }⩾0$ will rest against zero: it is profitable for us to decrease the sum of absolute values, that is, the sum $w_j^{+}+w_j^{-}$. More and more often it will turn out for some $j$ that both $w_j^{+}=0$ and $w_j^{-}=0$, but then $w_j=0$ too: the feature drops out of the model simply because its coefficient becomes zero. This is hard zeroing. Nothing of the kind happens in ridge regression: there the coefficients may tend to zero but never become exactly zero. This is the radical difference between $L_1$ and $L_2$ regularization.

Let us look at some regression problem: the coefficients $w_j$ as functions of the regularization coefficient, which here is also often called the selectivity parameter: the larger $\mu$, the stronger the regularization and the fewer features remain; selectivity strengthens the selection. In ridge regression, as the coefficient increases, all the $w_j$ begin to decrease, tend to zero, but never reach it. In LASSO, as $\mu$ increases, the coefficients of the features become zero one after another — that is how the noose works. In the end it has strangled all the coefficients; the model is over-regularized and has degenerated into a constant. At the start, with weak regularization, all the coefficients are nonzero, and the model is richest; in between are the intermediate variants. Like the regularization coefficient in ridge regression, $\mu$ in LASSO can be chosen, and there are methods for this — that is a separate story I am not telling. Some compromise variant will throw out part of the features; that is the main property of LASSO regularization. It is also called Tibshirani's lasso, after the author's surname. A good book where this can be found is The Elements of Statistical Learning, in effect the classic university textbook on machine learning, republished many times.

A geometric interpretation of feature selection

One more view of why $L_1$ regularization selects features. Take the equivalent statement: the least squares criterion is minimized under an inequality constraint on the $L_1$ or $L_2$ norm of the weight vector, which cannot exceed the threshold $\varkappa$.

The grey region is the region of admissible solutions. The ellipses are the level lines of $Q(w)$, and the global minimum of $Q$ is outside the allowed zone: it is that very unstable solution. So we look for the minimum only over the grey region the regularizer allows us. For $L_1$ it is a diamond, and look how interesting: wherever the quadratic form is located, the touching most often happens at a vertex of the diamond, on one of the axes — and there one of the coordinates becomes zero. For $L_2$ the region is a disc, and however we move the quadratic form around the space, only on a set of measure zero will we hit points where zeroing happens. In fact we will never hit them, among other things because of computational errors. Here is a vivid picture of why the $L_1$ norm leads to coefficients becoming zero, to feature selection.

One can generalize and take not necessarily $L_1$ or $L_2$ but other norms $L_p$, and look at what their equidistant surfaces look like — the points at a given distance from the origin. $L_0$ is simply the number of features used: we use one feature, either this one or that one. $L_{1/2}$ has a concave shape, and it will have greater selectivity: as the quadratic form is moved, touching at the vertices will happen more often. With $L_2$ and with $L_{\infty }$, the uniform norm, there will be no feature selection.

Non-smooth regularizers for selecting and grouping features

If the model is linear, the $L_1$, $L_2$ and even $L_p$ norms play the role of shrinkage. In general a regularizer is some function $R_{\mu }$, which one can try taking in different ways.

Elastic Net is a very useful thing, the sum of the $L_1$ and $L_2$ norms. Why is $\mu$ singled out separately here? Increasing it leads to feature selection, but among the selected features there may still be linearly dependent ones: we allow such multicollinearity but suppress it with $L_2$ regularization. The regularization coefficient and the selectivity parameter are slightly different parameters here, and the features are suppressed less aggressively. It can be done differently too: the methods I told you about at the end of the last lecture, the Support Feature Machine and the Relevance Feature Machine, define the regularization function differently. The main thing is that the regularizer must be non-smooth: if it has this sharp nose, however we construct it, there will be feature selection. Non-smoothness of the regularizer gives feature selection.

Principal component analysis

Let us devote the last third of the lecture to another problem with least squares, now — transforming the feature space. This is the third strategy for fighting multicollinearity. The first is to shrink the coefficient vector; the second, to zero some coefficients, throw out features; the third, to move to a smaller number of new features while preserving the information about the objects. This is a separate problem statement, and the answers $y$ are not needed at all in it: we work only with the «objects–features» matrix, and it does not matter to us whether we are going to solve a classification problem, a regression problem, or need a new feature representation of the objects for something else.

We want to find a space of smaller dimension that would represent the original objects more efficiently, with an explicit transition from the original space to the new one: the new features $g$ must be expressed explicitly through the x's. We want the reverse transition to be possible too: so that the old features can be reconstructed from the new ones, losing as little accuracy as possible. How to reconstruct — a model has to be set. If the model is linear, we get principal component analysis; with a nonlinear model we get other methods, which also exist, but we will not talk about them today. A linear combination of the new features with coefficients $u_{jt}$, unknown for now, should give, of course, not exactly the original features but their estimates ${\hat{f}}_j$, and we require the estimates to be as accurate as possible on all the objects of the training sample — least squares, a sum over all the features and all the objects. The criterion must be minimized simultaneously over the new features and over the coefficients $u_{jt}$, that is, we must find both the new features and the way to reconstruct the old ones from them.

The rows of $G$ multiplied by the columns of $U^{\mathsf{T}}$ — linear combinations of the new features — must give the old features: $G{U}^{\mathsf{T}}$ must be as close to $F$ as possible. Both matrices are unknown. We approximate $F$ by a product of two matrices of a smaller intermediate dimension $m$, so the rank of $G{U}^{\mathsf{T}}$ is at most $m$. Let us agree that we solve the problem to obtain $m$ linearly independent new features — then this is a matrix of rank $m<n$. If $F$ were also of rank $m$, the representation would be exact; but if all the matrices have full rank, it is a matter of approximating $F$ by a product of two matrices of lower rank, that is, of a low-rank matrix decomposition.

Look carefully at this set of conditions: you have already seen it somewhere today. The new features are also expressed linearly through the old ones, by the ready formula $G=FU$. If we write $U\Lambda =F^{\mathsf{T}}FU$ out column by column, we get ${\lambda }_t{u}_t=F^{\mathsf{T}}F{u}_t$ — the condition that $u_t$ is an eigenvector of the matrix $F^{\mathsf{T}}F$; the matrix notation says that the columns of $U$ are orthonormal eigenvectors of $F^{\mathsf{T}}F$. Similarly the columns of $G$ are orthogonal eigenvectors of the matrix $F{F}^{\mathsf{T}}$. Look how interesting: the columns of $G$ are the new features, and they are orthogonal as vectors of values on the objects of the training sample. That is why the transformation is also called orthogonalization.

And the last, important property is the accuracy of reconstruction. The squared norm of a matrix equals the sum of its eigenvalues, and the squared norm of the difference, the accuracy of the matrix decomposition, turns out to be the sum of the smallest eigenvalues. Now it is clear why the method is called principal component analysis: for the residual to be minimal and the accuracy of reconstruction maximal, we must take the largest eigenvalues, because the decomposition error equals the sum of the remaining ones. The largest eigenvalues and the eigenvectors corresponding to them are called the principal components: they contain the maximum of information about the original matrix $F$.

The link with the singular value decomposition

Up to notation, I think you have already seen the link with the singular value decomposition.

If the dimension is not reduced, the sum of the remaining ${\lambda }_j$ consists of zero terms, and the representation is exact. It coincides with the singular value decomposition if $G$ is represented as the product of an orthonormal matrix and a diagonal one with the roots of the eigenvalues on the diagonal. The transformation $U$ works both ways: $U^{\mathsf{T}}$ on the right of $G$ gives $F$, $U$ on the right of $F$ gives $G$, because $U^{\mathsf{T}}$ here is $U^{-1}$. And since the new features are uncorrelated, this is also the Karhunen–Loève transformation. Principal component analysis was probably reinvented dozens of times over the decades: it appeared in very different problems, in different notation, people derived everything anew and ended up arriving at it.

The effective dimension of a sample

When does principal component analysis bring benefit, and when is it useless? It brings benefit when there are linearly dependent features among the features. That means the sample is in fact distributed not over the whole space but lies in a subspace of smaller dimension. Principal component analysis moves explicitly into this subspace and finds a new basis in it — new features, the matrix $G$ — an orthogonal one even, which is very convenient. The original problem can now be reformulated in these features.

How can we tell whether the sample lies in a subspace of smaller dimension? With the same spectrum. The relative error $E_m$ — the squared norm of the difference divided by the squared norm of $F$, that is, the sum of the smallest eigenvalues divided by the sum of all — is a very convenient quantity from zero to one, and it can be plotted. If in one place a sharp drop has occurred, followed by a tail where the relative errors are already small and can be neglected, then the point of the drop $m$ is the dimension at which the error dropped sharply, and a further increase in the dimension of the subspace will not lead to a significant improvement in accuracy. So $m$ is the dimension of the space in which our sample is actually located.

Is it always so? Of course not. In many problems the plot is a monotonically decreasing function, decreasing to the end and never reaching values close to zero. Then the sample lives in a space of full dimension, and applying principal component analysis to it is useless and incorrect. We found the singular value decomposition, plotted the spectrum, and look at it with our eyes: is there a steep slope? If there is, then there are reasons why the sample ended up in a space of smaller dimension. For example, we constructed the features in such a way that some implicitly depend almost linearly on the rest, and we did not catch this dependence — we lacked the imagination to invent independent features. This happens fairly often, so it really does make sense to look at the dimension of the subspace the sample actually lives in.

Principal component analysis in linear regression

We have considered two problems: supervised learning, multiple linear regression, and unsupervised learning, principal component analysis. We solved both by least squares, understood their properties, and both are deeply tied to the singular value decomposition. Let us put a bold full stop here, take a deep breath — and now apply principal component analysis to solving the problem of multiple linear regression.

$U^{\mathsf{T}}w$ is the new parameter vector, now over the new features, and $Gz$ is our new model. We found $G$ and $U^{\mathsf{T}}$, threw the original $F$ away, and solve the least squares problem with the new feature matrix $G$; between the new and old coefficient vectors there is a linear link following from the properties of principal component analysis. And principal component analysis is expressed through the same singular value decomposition, only in it one has to take the $m$ largest eigenvalues and throw away the rest. The form of the solution is the same. Recall what ridge regression was for: to move the ${\lambda }_j$ away from zero, and we added $\tau$ to them. Now we move them away from zero differently: we simply throw the eigenvalues close to zero out of all the expressions as terms. This is equivalent to principal component analysis. An interesting link: both regularization and principal component analysis in the solution of the least squares problem are the same decompositions, in which we move the spectrum away from zero by different techniques.

The spectral method of least squares

This leads us to a mega-generalization of least squares. Having expressed least squares through the singular value decomposition, we understood how $L_2$ regularization is expressed through the SVD, and then that principal component analysis is expressed through the SVD. The structure of the solutions is absolutely identical; the only difference is how we move the small ${\lambda }_j$ away from zero. A generalization comes to mind at once.

We keep all the eigenvectors, we keep the structure of the solution. Ridge regression is the special case where $\tau$ is added to all the eigenvalues. But if an eigenvalue was large, why add $\tau$ to it? We need not: the main thing is that the small ones are moved away from zero. Principal component analysis is as if we added infinity to the small ${\lambda }_j$, that is, zeroed the corresponding terms in the sums where ${\lambda }_j$ is in the denominator. Something in between can also be done, in different ways. There are many variants of how to set the additions ${\delta }_j$; they can be invented and then the question posed: which of these spectral solutions suits us better for some reason? Such is the mega-generalization, not at all obvious: it would have been hard to guess it from the start. We arrived at it by comparing ridge regression with substituting principal component analysis into the solution of the least squares problem and noticing that the structure of these solutions is the same.

Low-rank decompositions: what comes next

And a spoiler of what lies ahead: in the second semester there will be a lecture on recommender systems, and a few more occasions to talk about low-rank decompositions. Principal component analysis is not the only statement.

Very often the matrix $F$ is given not in full but only on some subset of cells: missing values — for some objects not all features could be measured. The classical singular value decomposition assumes the matrix is fully known; there can be no gaps in it, and it is not even clear how to handle them. But there are more modern methods that solve such a problem explicitly: the sum of squares is still minimized, but only over the cells for which there are observations. It can be solved, for example, by the stochastic gradient method. This is of course no longer the SVD, and there will be no orthogonality properties for $G$ and $U$. The matrix $F$ is sparse, known only on a subset of values, and it is represented by a product of two dense matrices of lower rank, where there are no gaps any more. There are other variants too: a loss function that is not quadratic, additional non-negativity conditions. This is just a spoiler, so that you feel what other problem statements there are — things by no means end with principal component analysis.

Summary

Today we realised a very interesting role of the singular value decomposition. With it one can solve the problem of multiple linear regression; one can insert $L_2$ regularization, ridge regression, into it — also through the singular value decomposition; one can do principal component analysis — also through the singular value decomposition; and one can combine principal component analysis with the solution of the multiple linear regression problem and also get a solution through the singular value decomposition. Such an interesting and fairly general thing.

Seminar