# 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.

• @Pegah, you omited the dimensions of matrix $A$. Please check that I added them correctly. – mpiktas Jun 2 '11 at 14:38
• @Pegah, 2400 is not a large n. Such matrix will take about 40 Mb of RAM, which is not a lot for modern computers. – mpiktas Jun 2 '11 at 14:40
• @mpiktas thanks for the edits. I saw that recursive methods' time complexity is O(n!) in this link and Crout method is only O(n^3) – Pegah Jun 2 '11 at 14:43
• @Pegah $b = (X'X)^{-1} X'Y$ is just one way of saying that $b$ solves the equation $(X'X)b = X'Y$. Symbolically replacing $(X'X)\to A$, $X'Y \to B$, and $b \to x$ shows this is identical to $Ax = B$. – whuber Jun 2 '11 at 17:45
• @whuber thanks it is the answer. Could you please put it as an answer? So other users can see the answer. Thanks again. – Pegah Jun 2 '11 at 17:54

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$.

• hi, the problem is not with the size of n. the problem is that first I am not sure if my matrix in positive definite, so I can't use Cholesky decomposition. and second I don't know how to use Crout method because it needs an n by n matrix but my matrix is n by m. – Pegah Jun 2 '11 at 17:05
• @Pegah: I am confused. It is easy to prove that $X^TX$ is positive (semi-)definite for any $X$. Are you alluding to some numerical issues? – NPE Jun 2 '11 at 19:00
• On the other hand, forming $\mathbf X^T \mathbf X$ may not always be the best course of action (you square the condition number, among other things). At the very least, performing a QR decomposition of $\mathbf X$ might be better... – J. M. is not a statistician Aug 1 '11 at 7:01