Practical Linear Least Squares Calculations

In the previous installment, it was shown how a least squares problem can be set up to solve for linear model parameters. The specific example was to approximate a temperature versus resistance characteristic for an RTD sensor. The process involved:

• Collect independent "input" variables into a matrix.
• Collect the dependent "output" variable values into a column vector.
• Reformulate the overdetermined (and typically huge!) matrix problem as a (typically small!) Normal Equations problem.
• Solve the Normal Equations problem to obtain best-fit parameter values.

You should have no remaining fears about how that process works. But you might dread the prospect of some serious and tedious number crunching. Let us dispel that now.

 

Introducing Octave

What, Octave?^[1] Why not Matlab?^[2]

Go to the www.gnu.org Web site and download a copy of Octave that you can use on pretty much any kind of computer you want.

Now go to the www.mathworks.com Web site and try the same trick for Matlab. Be sure to verify that you aren't down $2000.00 US when you log out. ^[3] And don't forget to budget for the twice-annual upgrades.

It's not that Matlab isn't a great product. It's great in the sense of quality, great in the sense performance, great in the sense of comprehensive, great in the sense of authoritative, great in the sense of massive. Speaking of greatness, wouldn't it be great if those "great and benevolent ones" ^[4] granted me the entitlement to use such fabulous software without the price tag? (And sent me a workstation machine powerful enough to support it?) But for the rest of us, we don't need to bankrupt ourselves just to solve a few interesting math problems.^[5] So get a copy of Octave installed, and let's go calculating.

Oh, I should mention, the scripts for Octave and Matlab are so nearly compatible that anything you see here should work fine in either package.

 

Building an Octave script

Let's return to the example in the previous installment. The purpose was to fit a linear model with two parameters to calculate the temperature given a measurement of resistance for a nominally 10-Ohm copper RTD. The data was given in the form of a table. The first task is to enter the data into Octave. Use your favorite text editor to prepare a text file the following lines for the independent data values.

% Octave script 'rtdfit.m' to build linear model from copper RTD data.
% Set up "independent variable" (input) matrix with resistances and artificial constants.
MMat = [ 9.04     1;   9.23     1;   9.42    1;   9.65    1;  ...
         9.81     1;   10.00    1;  10.19    1;  10.39    1;  ...
        10.58     1;   10.77    1;  10.97    1;  11.16    1;  ...
        11.35     1;   11.58    1;  11.74    1;  11.93    1;  ...
        12.12     1;   12.32    1;  12.51    1;  12.70    1;  ...
        12.90     1 ]

The % notation indicates that the rest of the text line is comment. The [ and ] brackets indicate construction of a matrix. The matrix is entered by rows, with semicolon separators marking the end of each row. Use the "dot dot dot" notion to continue the data on the next text line. A final terminating semicolon after the closing bracket is omitted, which means display a copy of this item. After you establish correctness, you can add a trailing semicolon to omit unnecessary displays.

Now similarly set up the "dependent variable" vector.

% Set up "dependent variable vector" with corresponding temperatures.
VVec = [ 0  5  10  15  20  25  30  35  40  ...
         45  50  55  60  65  70  75  80  85  90  95  100 ]'

That final single-apostrophe character is no accident. I've done a little trick here... I set up the list of temperatures as a row vector (for Octave, a matrix with only one row) and then transposed it using a "tick" (single-apostrophe) operator to produce a column vector.

Now construct the normal equation matrices. To do this, use the matrix product MT · M. Matrix "multiplication" operations are expressed using the * operator, and you already know about transpositions using the "tick" operator.

% Set up the normal equation matrix terms
MTM = MMat' * MMat
MTV = MMat' * VVec
MTMinv = MTM^-1

That "caret minus one" notation specifies inversion of the matrix.

% Solve the normal equations to determine best fit parameters
VFit = MTMinv * MTV

Save this file as rtdfit.m and exit your text editor.

Now let's run the script to obtain the solution. Start the octave interpreter.

octave

In the octave environment, use the cd command to make the folder containing the script file your current working directory.

cd  'C:\projects\rtdfit'

Now type in the name of the script file you just set up, omitting the .m suffix, after the prompt.

>  rtdfit

Done. Amazing.

But let's not take anything on faith here. Let's apply this result to the original dataset to see what the model predictions are, and compare to the output data in the original data set. Add the following test to the end of your script file.

% Compute residuals to verify quality of fit
Modcurv = MMat * VFit
Fiterr = VVec - Modcurv

Now tell octave to run the script again. You will see the predictions of the model, and how close the model predictions are at each individual data point. To see the numerical results: look back to the previous installment.

 

More?

If you thought that was too easy, here is something that is even more amazing. The previous sequence was "by the book" and it took us a whopping half dozen commands to set up and solve the normal equations. But here is another way. Instead of the calculations shown previously, use the following.

 
% Alternate LLS parameter computation by pseudo-inverse method
VFit = MTM \ VVec

In effect, Octave "knows" what you have in mind, and constructs for you the equivalents of the Normal Equation calculations automatically. It is almost too magical.

 

Easy extensions

Let's consider some additional twists that commonly occur.

• We covered the case of a constant and a single input. What if there are more input variables? Just add an additional matrix column for each additional input variable, and there will be one new parameter variable associated with it in the solution vector.

• What happens if input variables multiply each other? The LLS method really doesn't care how input vectors were constructed. Set up a new matrix column containing the values of the products of the two variables, and there will be one new parameter variable associated with the product terms in the solution vector.

• What happens if there are two output variables that the model is supposed to predict, not just one? No problem. Instead of a single column vector, expand the VVec term to be two vectors. The best fit will become a matrix with two columns of parameter values, the first column related to producing the first output, the second column related to producing the second output.

Before we leave the topic of Linear Least Squares, there are a few more things you probably need to know. However, before diving too deeply into secondary considerations, let's experiment with building Kalman filter dynamic equations using what we have so far.

 

---