CODEWITHSUNDEEP
pythonmatlablinear-regressionlarge-data
PeterC
PeterC

Reputation: 21

Generalized least square on large dataset

I'd like to linearly fit the data that were NOT sampled independently. I came across generalized least square method:

b=(X'*V^(-1)*X)^(-1)*X'*V^(-1)*Y

The equation is Matlab format; X and Y are coordinates of the data points, and V is a "variance matrix".

The problem is that due to its size (1000 rows and columns), the V matrix becomes singular, thus un-invertable. Any suggestions for how to get around this problem? Maybe using a way of solving generalized linear regression problem other than GLS? The tools that I have available and am (slightly) familiar with are Numpy/Scipy, R, and Matlab.

Upvotes: 2

Views: 2304

Answers (2)

Pursuit
Pursuit

Reputation: 12345

Instead of:

b=(X'*V^(-1)*X)^(-1)*X'*V^(-1)*Y

Use

b= (X'/V *X)\X'/V*Y

That is, replace all instances of X*(Y^-1) with X/Y. Matlab will skip calculating the inverse (which is hard, and error prone) and compute the divide directly.

Edit: Even with the best matrix manipulation, some operations are not possible (for example leading to errors like you describe).

An example of that which may be relevant to your problem is if try to solve least squares problem under the constraint the multiple measurements are perfectly, 100% correlated. Except in rare, degenerate cases this cannot be accomplished, either in math or physically. You need some independence in the measurements to account for measurement noise or modeling errors. For example, if you have two measurements, each with a variance of 1, and perfectly correlated, then your V matrix would look like this:

V = [1   1; ...
     1   1];

And you would never be able to fit to the data. (This generally means you need to reformulate your basis functions, but that's a longer essay.)

However, if you adjust your measurement variance to allow for some small amount of independence between the measurements, then it would work without a problem. For example, 95% correlated measurements would look like this

V = [1      0.95; ...
     0.95   1  ];

Upvotes: 2

duffymo
duffymo

Reputation: 308763

You can use singular value decomposition as your solver. It'll do the best that can be done.

I usually think about least squares another way. You can read my thoughts here:

http://www.scribd.com/doc/21983425/Least-Squares-Fit

See if that works better for you.

I don't understand how the size is an issue. If you have N (x, y) pairs you still only have to solve for (M+1) coefficients in an M-order polynomial: 

y = a0 + a1*x + a2*x^2 + ... + am*x^m

Upvotes: 1

Related Questions