Creating simple procedure for The Least-Square $m^\text{th}$ Degree Polynomials

Question

I am CS major, taking Computational Numerical Analysis course. Instructor gave us freedom of choice, we were allowed to use anything or any computer language we picked, I picked Mathematica.This is my 3rd day with the program.

I am writing a simple procedure to calculate the Least-Square $m^\text{th}$ Degree Polynomials based on data sample. Algorithm is somewhat simple: http://www.efunda.com/math/leastsquares/lstsqrmdcurve.cfm
This is how far I got:

y_sample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
x_sample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}(*Some sample data*)

f[x_] = Sum[Subscript[a, i] x^i, {i, 0, 4}](*generating polynomial, in this case, degree of 4*)    
Sum[(y_sample[[i]] - f[x_sample[[i]]])^2, {i, 0, 4}] 

Last line of code should roughly correspond to

$$ \Pi = \sum^n_{i = 1}\left[y_i - f(x_i)\right]^2 = \sum^n_{i=1}\left[y_i- (a_0 + a_1 x_i + a_2 x_i^2 + \ldots + a_m x_i^m)\right]^2 = \text{min} $$

List of problems and questions, I am facing:

1. Executing code above, produces number of errors that I don't mean much to me.
2. How do I differentiate in respect to all $\mathbb{a}_n$
3. Differentiation should produce number of linear equations, how do I set the results of differentiation to zero and solve for unknown coefficients.
4. Outputs are enormously long, is there a way to suppress some outputs !

Answer 1

You have some errors in your syntax:

• you name your lists x_sample and y_sample, but in Mathematica, an underscore is not allowed in names (as it is reserved to patterns).
• your last sum runs from 0, but in Mathematica, the first element in a List has index 1
• your last sum should run until the number of data points, not 4
• furthermore, I would advise you to use SetDelayed(:=) in your function definition, as it prevents possible naming conflicts
• also in your function definition, letting the sum run from 0 will give problems for x=0, as it produces 0^0, which is Indeterminate.

This gives:

ysample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
xsample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}

f[x_] := Subscript[a,0] + Sum[Subscript[a, i] x^i, {i, 1, 4}]
result = Sum[(ysample[[i]] - f[xsample[[i]]])^2, {i, 1, Length[ysample]}]

Next, to calculate the derivative with respect to, say, Subscript[a,0], use D[result,Subscript[a,0]].

Last, to get a solution to the set of linear equations, use Solve:

sols = Solve[Table[D[result, Subscript[a, i]] == 0, {i, 0, 4}], Table[Subscript[a, i], {i, 0, 4}]]

To put the values of the coefficients in your polynomial, use Replace(/.):

final[x_] := (f[x] /. sols)[[1]]
final[x]
final[2]

EDIT

As you can see in the result, there are some values which are almost zero (10^-11 or smaller). These are simply the result of some internal round-off errors, and should of course be equal to zero. You can easily add a piece of code to discard all coefficients which are a factor smaller than the largest coefficient.

The largest coefficient is given by:

M = Max[Apply[#2 &, sols, {2}][[1]]

(don't mind the syntax too much yet, it is just a way to select the values from sols in order to be able to use Maxon it)

Next we make the replacement for every x which is at least 10^6 times smaller than the maximum:

sols = sols/. x_ /; x < M / 10^6 :> 0
final[x_] := (f[x] /. sols)[[1]]
final[x]
final[2]

EDIT 2

It appears Mathematica has a built-in function to deal with near-zero numerics: Chop. It's al lot more easy to use than the method I proposed in my former edit; just write Chop[sols] and it will automatically thread over the rules. Thanks J.M. for pointing that out!

Then the final code would be:

ysample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
xsample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}

f[x_] := Subscript[a,0] + Sum[Subscript[a, i] x^i, {i, 1, 4}]
result = Sum[(ysample[[i]] - f[xsample[[i]]])^2, {i, 1, Length[ysample]}]

sols = Solve[Table[D[result, Subscript[a, i]] == 0, {i, 0, 4}], Table[Subscript[a, i], {i, 0, 4}]]

final[x_] := (f[x] /. Chop[sols])[[1]]
final[x]
final[2]

Answer 2

freddieknets already told you what you're supposed to be doing. Here's how I'd have done your approach to least squares:

n = 4;
f[x_] = Map[C, Range[0, n]].x^Range[0, n];


Solve[Thread[
    D[Total[(ysample - f /@ xsample)^2], {Map[C, Range[0, n]]}] == 0],
    Map[C, Range[0, n]]] // First // Chop
{C[0] -> 0, C[1] -> 10., C[2] -> 0, C[3] -> 0, C[4] -> 0}


f[x] /. %
10. x

Notes:

• C was always intended for indexed "arbitrary constant" use; this is one case where it certainly is useful.

• The dot product is useful for assembling a polynomial from its coefficients. An alternative is to use a "Horner approach": f[x_] = Expand[Fold[(#1 x + #2) &, 0, Map[C, Range[n, 0, -1]]]]

• The least-squares objective function is compactly implemented here as Total[(ysample - f /@ xsample)^2]; Sum[] is useful notationally, but in Mathematica, there are always other methods.

• Remember that in an optimization such as this, one is supposed to be solving for the zeroes of the gradient; in this case, note that D[f[x, y, (* other variables *)], {{x, y, (* other variables *)}}] is the way to quickly generate a gradient.

• Thread[list == 0] equates all the components of list to 0.

• As I told freddieknets, Chop[] is very useful for problems such as this, when you want to treat tiny terms as actually being zero.

---

For giggles, here's one way I'd have done, if none of Fit[], FindFit[], LeastSquares[], or PseudoInverse[] were available:

n = 4;
vm = Take[LinearAlgebra`VandermondeMatrix[xsample], n + 1];
{qm, rm} = QRDecomposition[Transpose[vm]];
LinearSolve[rm, qm.ysample].x^Range[0, n] // Chop

I could have used CholeskyDecomposition[] (the normal equations route), but QRDecomposition[] is often more numerically sound for applications like these.

There's also SingularValueDecomposition[], but that requires slightly more machinery to discuss properly.

Finally, I gave a rundown of most of the usual methods for linear least squares in this stats.SE answer; you might be interested in it.

Answer 3

Also, remember you can check your results with:

ysample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
xsample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}(*Some sample data*)
model = a x^3 + b x^2 + c x + d;
fit = FindFit[data = Transpose@{xsample, ysample}, model, {a, b, c, d}, x]
modelf = Function[{x}, Evaluate[model /. fit]]
Plot[modelf[x], {x, 0, 1}, Epilog -> Map[Point, data]]