profile for Kells1986 at Stack Overflow, Q&A for professional and enthusiast programmers

Friday, 9 March 2012

Moore-Penrose Pseudo-inverse Best Fit

Football Fan? Why not check out www.playerrank.co.uk to vote for your favourite ever player?

Before we get started I just want to say that I've signed this blog up to technorati, in order to get some more hits. I need to "claim" this blog by posting the following code: HMR6GJ2P49B6

Now that's done let's get to the meat of today's blog...

I was in the pub on Monday with a few of my colleagues, where we were discussing our day's work. Two colleagues had been doing least-squares best fitting and one was complaining that the code she used wasn't working.

She was trying to fit an exponential curve to a data set using some kind of Chi-Squared minimizing function in IDL and it was all going horribly wrong. At this point I asked why she wasn't using the Moore-Penrose Pseudoinverse matrix to do a linear fit on the logged data values. To my surprise none of my colleagues had heard of this technique, so I thought it would make a good blog.

Suppose you had the following system of equations:

There is no one point where the lines all intercept, since the three equations are linearly independent, so a least squares solution is sought.

Decompose this system into matrix form, Ax=b. You now have a system that looks like this:

A = [2  -1]  x = [x] b = [ 2]
      [1   1]        [y]       [ 5]
      [6  -1]                   [-5]

To find the least squares value of x and y all you need to do is the following:


in this case the result is  x = -0.094595 and y = 2.445946. 

In IDL the code to do this is really easy:

x = (invert(transpose(A)##A)##transpose(A))##B


Anyway, that's all for today's blog, if you've got any questions please feel free to comment below, it would be good to hear from you.


  

4 comments:

  1. Thanks, this has been very helpful. I'm jumping in the deep end, and I've seen a lot of code, and a lot of theory equations, but this linked them up nicely. Well done!

    ReplyDelete
  2. Hey, thanks. I've kind of forgotten about about this blog but it seems to be getting hits. I'm gonna take some time over the next couple of weeks to put up some more content. I'll also try and respond if anyone has any specific maths/coding questions.

    ReplyDelete
  3. That's a really nice explanation. Could you recommend some c/c++ or C# code that is able to achieve this using matrix operations.

    I know it's taking the easy option but I have a very tight deadline and if there is an example out there I would be extremely grateful. I would also ensure that I understood the code before I released it to the world.

    ReplyDelete
  4. I don't have any code written myself to do this in C, but I would suggest checking Numerical Recipes, or having a look at the LAPACK library (http://www.netlib.org/lapack/index.html#_standard_c_language_apis_for_lapack)

    ReplyDelete