# Least Square Data Fit /Design Matrix

I am trying to learn how to use the DesignMatrix function in Mathematica and how to perform a LeastSquares fit. I am at end of my rope and seek expert intervention.

Say, I have some experimental data and some experimental fit coefficients

size = 20; (*no of points *)
expdata = (Array[#, size] & /@ {x, y, z})~Join~{ConstantArray[1, size]};
coeffs = {{a1, b1, c1, d1}, {e1, f1, g1, h1}, {i1, j1, k1, l1}}; (*can be larger*)


My experimental equations would look like:

lhs = coeffs.expdata; (*Left hand side *)
rhs = Array[#, size] & /@ {u, v, w}; (*Right hand side *)
eqs = Flatten[ Table[ lhs[[#1, i]] == rhs[[#1, i]] & /@ Range, {i, 1, size}]];


I know I can solve these equations exactly using the Solve function as follows for a small number of data points:

Solve[ eqs, Flatten[coeffs]]


But assuming I have a large number of noisy data points, how do I setup a LeastSquares problem. In specific, I’d like to understand how to use the DesignMatrix function….

• You've three disjoint problems in there, one for each sublist in coeffs – Dr. belisarius Jan 28 '14 at 20:30

Let's say you have a set of data:

data = {{-2., -9.424}, {-1.5, -2.586}, {-1., -3.047}, {-0.5, -1.203}, \
{0., 0.551}, {0.5, 4.566}, {1., 12.077}, {1.5, 21.2118}, {2.,
44.752}};


and you suspect it can be modeled with a 3rd order polynomial. You can set up the design matrix using the following:

a = DesignMatrix[data, {x, x^2, x^3}, x]


This link provides some useful information on how the least squares approximation works. Simply, if the equation a x = b has no solution, multiply both sides by the Transpose of a and solve the resulting equation. This can be done in Mathematica by:

soln = Solve[Transpose@a.a.{k1, k2, k3, k4} == Transpose@a.data[[All, 2]],
{k1,k2, k3, k4}]
(* {{k1 -> 0.368633, k2 -> 3.98313, k3 -> 4.23867, k4 -> 2.3121}} *)


Which gives us a decent fit

Plot[k1 + k2 x + k3 x^2 + k4 x^3 /. soln, {x, -2, 2},
Epilog -> {Red, PointSize[0.02], Point /@ data}] which is also identical to doing it the simple way

nlm = NonlinearModelFit[data,
k1 + k2 x + k3 x^2 + k4 x^3, {k1, k2, k3, k4}, x]
nlm["BestFitParameters"]
(* {k1 -> 0.368633, k2 -> 3.98313, k3 -> 4.23867, k4 -> 2.3121} *)


Note I use NonlinearModelFit instead of LinearModelFit here only because the output more closely resembles the output of soln

• What a nice, complete answer! +1 – ciao Jan 29 '14 at 0:37
• Thank you. This is a much better explanation that the Help pages on Ma. – Pam Jan 29 '14 at 2:01