4

This is the equation I am trying to solve:

h = (X'*X)^-1*X'*y

where X is a matrix and y is a vector ((X'X)^-1 is the inverse of X-transpose times X). I have coded this in Matlab as:

h = (X'*X)\X'*y

which I believe is correct. The problem is that X is around 10000x10000, and trying to calculate that inverse is crashing Matlab on even the most powerful computer I can find (16 cores, 24GB RAM). Is there any way to split this up, or a library designed for doing such large inversions?

Thank you.

Jordan
  • 3,998
  • 9
  • 45
  • 81
  • How large of a matrix? Can you show me the vector y and matrix X. Or at least give me an idea of what it is like? – Derek W Mar 12 '13 at 23:49
  • What is the dimension of `X` and `y`? You might want to compute an approximation of the inverse if the dimension is very high. Don't you miss some `*`s btw? – sfotiadis Mar 12 '13 at 23:50
  • Edited question with the info. X is about 10,000x10,000. The dimension of the vector y matches (also about 10,000 long). – Jordan Mar 12 '13 at 23:51
  • do you know anything about the structure of the matrix? like if its sparse etc ? BTW I reckon instead of '\' you can use inv() function as '\' solves the least square stuff (pseudo inverse) – amas Mar 12 '13 at 23:53
  • Is is sparse? do you need to keep the `double` precision? do you have access to a gpu? – bla Mar 12 '13 at 23:54
  • X is a special strucure, it's kind of like 3000-diagonal (first column will be 3000 non-zero then 4000 zero, second column will be 1 zero then 3000 non-zero then 3999 zero, etc. To correct my earlier comment, X is actually 7,000x10,000, so X'*X is 10,000x10,000. – Jordan Mar 12 '13 at 23:55
  • This will be probably useful to you... http://math.stackexchange.com/questions/16940/in-place-inversion-of-large-matrices – bla Mar 12 '13 at 23:56

4 Answers4

1

That looks like a pseudo inverse. Are you perhaps looking for just

h = X \ y;
Ben Voigt
  • 277,958
  • 43
  • 419
  • 720
0

I generated a random 10,000 by 10,000 matrix X and a random 10,000 by 1 vector y.

I just broke up my computation step by step. (Code shown below)

  1. Computed the transpose and held it in matrix K
  2. Then I computed Matrix A by multiplying K by X
  3. Computed vector b by multiplying K by vector y
  4. Lastly, I used the backslash operator on A and b to solve

I didn't have a problem with the computation. It took a while, but breaking up the operations into the smallest groups possible helped to prevent the computer from being overwhelmed. However, it could be the composition of the matrix that you are using (ie. Sparse, decimals, etc.).

X = randi(2000, [10000, 10000]);
y = randi(2000, 10000, 1);
K = X';
A = K*X;
b = K*y;
S = A\b;
Derek W
  • 9,708
  • 5
  • 58
  • 67
  • 3
    randi is creating a single value in the range 1-10000, not a matrix 10,000x10,000. The final solution should be a vector. – Jordan Mar 13 '13 at 00:13
  • My mistake! Totally looked up `randint` syntax instead! I'll update my answer. – Derek W Mar 13 '13 at 00:17
  • Thanks. There is a small issue with the way I need to implement the code though: I need to apply this same X to many different values for the y vector, and I don't want to have to the invert each time. So, y can't be introduced until the very last step. – Jordan Mar 13 '13 at 01:09
0

If you have multiple machines at your disposal, and you can recast your problem into the form h = X\y as proposed by @Ben, then you could use distributed arrays. This demo shows how you can do that.

Edric
  • 23,676
  • 2
  • 38
  • 40
0

Jordan, Your equation is exactly the definition for "Moore-Penrose Matrix Inverse".

Check: http://mathworld.wolfram.com/Moore-PenroseMatrixInverse.html

Directly using h = X \ y; should help. Or check Matlab pinv(X)*y

Dazhao Liu
  • 1
  • 1