I'm having some troubles while fitting data with high order polynomials (typically order 6). I looked at the polyfit function in 'matplotlib/mlab.py' and I must say I don't understand why the least-square problem isn't solved using the matrix multiplication instead of the * multiplication. Unless I misunderstood something, shouldn't the following line:
c = array(linear_algebra.inverse(Xt*X)*Xt*y) # convert back to array
be rewritten like this?
c = array(matrixmultiply(matrixmultiply( linear_algebra.inverse(matrixmultiply(Xt,X)), Xt), y) ) # convert back to array