# Practical Linear Least Squares Calculations

Developing Models for Kalman Filters

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 ` M`

. Matrix "multiplication"
operations are expressed using the ^{T} · M`*`

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.

Footnotes:

[1] *Octave* at www.gnu.org

[2] *Matlab* at www.mathworks.com

[3] There are steeply discounted academic packages. They remind me
of the lovely ballad by Tom Lehrer...

He gives the kids free samples

because he knows full well

that today's experimenters

will be tomorrow's clientÃ¨le...

(*The Old Dope Peddler*, from *Songs by Tom Lehrer*,
Reprise 6216.)

[4] As Mitt Romney so clearly explained to us,
*corporations* are people too. Except not the benevolent part —
that's illegal. Look up *Dodge vs. Ford Motor Company, 1919.*

[5] Free as in *freedom* is not the same as
free as in *beer.* Support the software you use at
Octave - Donate.