{"id":81,"date":"2013-09-23T00:56:59","date_gmt":"2013-09-23T06:56:59","guid":{"rendered":"http:\/\/inductivebias.com\/Blog\/?p=81"},"modified":"2014-05-05T04:03:30","modified_gmt":"2014-05-05T10:03:30","slug":"least-squares-curve-fitting","status":"publish","type":"post","link":"http:\/\/inductivebias.com\/Blog\/least-squares-curve-fitting\/","title":{"rendered":"Least Squares Curve Fitting"},"content":{"rendered":"
\"LeastSquaresPoly5\"<\/a>Today we’re talking about linear least squares curve fitting, but don’t be fooled by the name. I’ll show how to use it to fit polynomials and other functions, how to derive it, and how to calculate it efficiently using a Cholesky matrix decomposition.<\/strong><\/span><\/p>\n

\nWe’ve got one set of input vectors where each one also has a corresponding output vector, and we’ve got a second set of input vectors that don’t yet have outputs. We want to learn the relationship between inputs and outputs using the first group and then use it to predict the outputs for the second group. This is your classic supervised learning problem, but unlike the Naive Bayes algorithm in the last post, this time we want to get actual numbers, not just which class the input belongs to.\u00a0 We can do this by fitting a line or a polynomial to our data, and both of these functions can be fit using linear least squares.\n<\/p>\n

Derivation<\/span><\/h2>\n

Let’s say we have an input vector X, and we want to find the linear function that best predicts an output scalar y. That is y = C*X (C dot X) where C is the vector of coefficients we’re solving for. In reality we have several X’s each with it’s own y, so X becomes a matrix with each row being a single input vector and Y becomes a vector of all of the outputs for each input vector.<\/p>\n

\\begin{bmatrix} X_1 \\\\ X_2 \\\\ ... \\end{bmatrix} C = \\begin{bmatrix} y_1 \\\\ y_2 \\\\ ... \\end{bmatrix}<\/p>\n

If the number of training points is exactly the same as the length of an input vector of a training point then the X matrix will be square and all we have to do is solve the resulting system of linear equations for the C vector.\u00a0 Barring some sort of matrix singularity, this C vector will produce a function which goes perfectly through all of the training points. Realistically, the number of training points (think thousands) will far exceed the number of input variables (think dozens). It’s usually not possible to pick a C which perfectly hits every training point, so we have to pick one that fits them as well as possible. In the case of least squares this means picking the one that gets the least sum of squares of the error for each training point, hence the name of the algorithm.<\/p>\n

\\underset{C}{\\arg\\min} (XC-Y)^2<\/p>\n

That is, we want to find the C vector that gives the minimum value if we multiplied it by our X matrix and subtracted the desired output and squared it. How do we do this? It’s simple calculus. Remember Fermat’s theorem<\/a>? Probably not by name, but it’s the thing where you set the derivative equal to zero to find minimums and maximums.<\/p>\n

\u00a0f(X) = (XC-Y)^2\\\\ {{\\partial f} \\over {\\partial C} } = 2{X^T}(XC-Y) = 0\\\\ 2{X^T}XC-2{X^T}Y = 0\\\\ {X^T}XC = {X^T}Y <\/p>\n

This looks exactly like our previous formula except with an X transposed prefacing each side. If X is an n by m matrix where n is the number of data points and m is the length of an input vector then X^T * X will be m by m and X^T * Y will be m by 1. This is good because it means the system of equations we need to solve to find C is square, so it has a unique solution, but its size is also independent of the number of training\u00a0 points. In addition any matrix of the form X^T * X has a few other useful properties we can take advantage of. Specifically, it is symmetric (that is Mij = Mji), and it is positive definite (that is a^T X a > 0 for all a). The savings from symmetricness is fairly straightforward: almost half of the matrix doesn’t need to be calculated or considered because it’s going to be the same as the other half. The positive definiteness guarantees that out solution is a minimum, and it turns out to be a requirement for the method we’ll use.<\/p>\n

Solving the System of Equations
\n<\/span><\/h2>\n

At this point we have a system of linear equations which we could solve using your generic linear algebra class Gaussian elimination. However, that’s not a very good idea for two reasons: it’s\u00a0unstable and it’s slow.\u00a0 Unstable in this case meaning the loss of precision from using computer arithmetic will be compounded by the operation, and slow meaning, well, slow.<\/p>\n

There are two algorithms I consider good candidates for solving this problem: Cholesky decomposition and QR decomposition.\u00a0 Both of these techniques work by breaking the matrix into two matrices that produce a problem which is easier (i.e. computationally cheaper) to solve. The primary difference being that the Cholesky runs on the already calculated X^T*X, while the QR computes its factorization of X^T*X directly from X itself.\u00a0 The QR is a bit more stable since some accuracy is lost when computing X^T*X, but the Cholesky is faster and in my experience plenty stable enough. We’ll be using the Cholesky today. Let our X^T * X be A, and the Cholesky decomposition is arranged as follows.<\/p>\n

\"WikipediaCholesky\"<\/a><\/p>\n

The important thing to note here is that L is lower diagonal. This means solving a system of equations with it is as simple as back-substitution.\u00a0 I’m not going to go much into the derivation of the decomposition itself, but suffice to say you could set the multiplied out form equal to A and work out the formulas. One of which looks like the following:<\/p>\n

\"WikipediaCholeskyFormula\"<\/a><\/p>\n

If you looked at this and immediately thought “well, that seems silly to take the square root and then just end up multiplying it back together in the next step” then congratulations you’re smarter than I am. I’m pretty sure I had to be told that it’s possible to factor out the diagonal terms in a separate matrix D and avoid those square-roots entirely without really adding any extra work. This form is called the LDL decomposition, and it looks something like this:<\/p>\n

\"LDL\"<\/a><\/p>\n

\"LDLformula1\"<\/p>\n

\"LDLFormula2\"<\/a><\/p>\n

This is what I recommend using to solve linear least squares problems. Calculating the decomposition is as simple as just plugging these formulas into your code, but you do need to be a bit careful as both formulas make use of earlier calculations of each other. You have to walk through and calculate the D’s for each row at the same time you calculate a row of L.<\/p>\n