Solving a regression problem I wanted to solve such a regression problem: 
$$Y = Xb + e$$
where $X$ is a $m$ by $n$ matrix, resulting in: b = (X'X)-1X'Y as a solution. Since $n$ is quite large (2400), I can't use the conventional methods to calculate the inverse of $X'X$. So I thought about using LU Decomposition using Crout method. In this link an implementation of this method is available which requires an $n$ by $n$ matrix $A$ to solve the equation $Ax=B$. Does it mean that I should use $X'X$ matrix as the input ($A$)? Please note that the problem is that it requires an $n$ by $n$ matrix instead of an $m$ by $n$ matrix. 
 A: You say that you need to solve an ordinary least squares problem on 2400 variables.
There are two assumptions that I think you need to revisit:
Assumption 1: that you need to compute the inverse of $X^TX$.
Assumption 2: that solving ordinary least squares on 2400 variables requires specialized methods.
I'll examine them in turn:
Assumption 1: that you need to compute the inverse of $X^TX$.
A better way to solve OLS using normal equations is by computing the Cholesky factorization of $X^TX$. See section 5.3.2 of Golub and van Loan for details. They state that the entire algorithm, including computing $X^TX$, $X^Ty$ as well as performing the Cholesky factorization and back-substitution, requires $(m+n/3)n^2$ floating-point operations.
Assumption 2: that solving ordinary least squares on 2400 variables requires specialized methods.
You don't say what kind of hardware you have at your disposal, so I'll assume that you have access to a typical mainstream PC.
First of all, a 2400x2400 of 64-bit floats requires just 44MB of memory.
Secondly, computing Cholesky decomposition of a matrix of this size takes half a second on my desktop PC using Numerical Python (numpy). This is the dominant computation once you have computed $X^TX$ and $X^Ty$.
