7
$\begingroup$

Is there a way to speed up $Z(Z'Z)^{-1}Z'$ type matrices? I am implementing the expression below directly using a matrix language and my program frequently crashes while if I run OLS on them using a pre-fabricated command, it is not an issue.

Is there a tip you guys might have to compute these matrices efficiently?

The goal here (but that is just an aside) is to implement the following estimator \begin{eqnarray} (X' P X - \sum_{i=1}^{n} P_{ii} X_{i}X_{i}' - \alpha X'X)^{-1} (X' P_Z y - \sum_{i=1}^{n} P_{ii} X_{i} y_{i} - \alpha X' y) \end{eqnarray} Now, $\alpha$ is the smallest eigenvalue of $(\overline{X}'\overline{X})^{-1} (\overline{X}' P_Z \overline{X} - \sum_{i=1}^{n} P_{ii} \overline{X}_{i} \overline{X}'_{i})$ where $\overline{X} = [y,X]$.

I am pretty confident that once I have found an efficient way to compute the projections, I can easily implement the rest easily.

Thanks so much!

Improve this question
$\endgroup$
10
  • 4
    $\begingroup$ The matrix $X'X$ is symmetric which can be used to speed up the inverse computation (e.g. compute Cholesky decomposition, then use it to compute the inverse). There is not much to do to speed up the matrix product ; first compute $X(X'X)^{-1}$, then $X(X'X)^{-1}X'$ (you know that it is symmetric so you can compute only e.g. the lower triangular part). $\endgroup$
    – Elvis
    Commented Mar 2, 2015 at 12:34
  • $\begingroup$ what do you mean 'implementing them verbatim'? $\endgroup$
    – user603
    Commented Mar 2, 2015 at 12:41
  • $\begingroup$ "from scratch" I guess :D $\endgroup$
    – Elvis
    Commented Mar 2, 2015 at 12:48
  • 3
    $\begingroup$ Do you need the entire hat matrix for anything in particular? Or are you only after part of it (such as the diagonal)? The regression program you use may already calculate the QR decomposition, which may be useful if you can get it. $\endgroup$
    – Glen_b
    Commented Mar 2, 2015 at 13:07
  • 1
    $\begingroup$ Hirek, I think it makes more sense to take your question back to where it was when it got its two answers, which deal fairly well with a perfectly good question, and then to ask a new question which reflects more nearly your actual problem. [It's not necessary to make any actual edits to your question to achieve that, you can simply roll it back in the edit history.] $\endgroup$
    – Glen_b
    Commented Mar 2, 2015 at 21:33

2 Answers 2

Reset to default
6
$\begingroup$

Using QR decomposition (which ought to be available if you already have calculated the regression):

Let $X$ have $n$ rows and $p$ columns and be of full column rank.

$H=X(X'X)^{-1}X'$

$=QR(R'Q'QR)^{-1}R'Q'$

$=QR(R'R)^{-1}R'Q'$

But if $R_1$ is the first $p$ rows of $R$ then $R'R=R_1'R_1$

$=QR(R_1'R_1)^{-1}R'Q'$

Now let $Q=(Q_1,Q_2)$ where $Q_1$ is the first $p$ columns of $Q$. Then $QR=Q_1R_1$.

$=Q_1R_1R_1^{-1}(R_1')^{-1}R_1'Q_1'$

$=Q_1Q_1'$

Where $Q_1$ is $n\times p$.

So if you have the QR decomposition of $X$, then the hat matrix is fairly simple.

Note that some regression programs will give $Q_1$ automatically. [It's also possibly that a regression program will have performed pivoting. That shouldn't impact the calculation of the hat matrix though.]

Improve this answer
$\endgroup$
23
  • $\begingroup$ Does QR work for non-square matrices as well? @Glen_b $\endgroup$
    – Hirek
    Commented Mar 2, 2015 at 13:40
  • $\begingroup$ What do you mean "work" -- which matrix is non-square? In regression, when we take $X=QR$, we have $X$ is nxp, $Q$ is nxn, and $R$ is nxp but upper triangular (so only the top pxp part is nonzero). $\endgroup$
    – Glen_b
    Commented Mar 2, 2015 at 13:40
  • $\begingroup$ X is n by k (as is usual in regression analysis) so X'X is a square. $\endgroup$
    – Hirek
    Commented Mar 2, 2015 at 13:42
  • 2
    $\begingroup$ @Elvis Actually, I believe the QR itself is more expensive than Cholesky ... but the subsequent calculation of the hat matrix is cheap, and it's considerably more widely used for regression since it's more stable. The question as it stood when I answered it implied the problem was being run after a regression. $\endgroup$
    – Glen_b
    Commented Mar 2, 2015 at 20:28
  • 2
    $\begingroup$ @Hirek What kind of pivoting is being used? Column pivoting? Or something else? [Having already changed from a simple regression question about a hat-matrix to something quite different, it now seems to be moving from a statistical problem to a software problem (I don't understand how to make a particular piece of software do something). It should probably roll back to where question and answers corresponded.] $\endgroup$
    – Glen_b
    Commented Mar 2, 2015 at 20:41
3
$\begingroup$

Try using the SVD. E.g., since $$ X = U\Sigma V^T, $$ then $$ (X^T X)^{-1} = (V\Sigma^2 V^T)^{-1} = V\Sigma^{-2} V^T, $$ and thus $$ X(X^T X)^{-1}X^T = U I_r U^T = U_r U_r^T, $$ where $I_r$ is an $n\times n$ identity matrix with $r\leq n$ ones on the diagonal (upper part), and $n-r$ zeros on the lower diagonal, where $r$ is the rank of $X$.

This will likely speed up your computation.

Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ I doubt it.The SVD is not a cheap computation. $\endgroup$
    – Elvis
    Commented Mar 2, 2015 at 12:48
  • $\begingroup$ Thanks @TommyL but I cannot quite follow your derivation. Is the last one an equality of the entire projection matrix, right? $\endgroup$
    – Hirek
    Commented Mar 2, 2015 at 13:07
  • 2
    $\begingroup$ Yes, that's correct. I've updated the answer to make this clear. @Elvis is right, though, the SVD is something like $\mathcal{O}(n^3)$, so it is not cheap, but if you repeatedly need to project, you may gain a lot. In particular if $n<<p$ or $p<<n$, you can save greatly. Note also that $UI_r U^T$ corresponds to selecting the $r$ first columns of $U$ and compute $U_rU_r^T$. $\endgroup$
    – Tommy L
    Commented Mar 2, 2015 at 13:26
  • 3
    $\begingroup$ @Tommy I think QR is also $O(n^3)$; I don't think any stable algorithm for computation of $X(X'X)X'$ will be of faster leading term than that; the only difference may be the factor out the front. Given the way the question is written now, this might turn out to be a better choice. $\endgroup$
    – Glen_b
    Commented Mar 2, 2015 at 21:14
  • 1
    $\begingroup$ You’r right, the asymptotic complexity is the same for all these methods. It must take into account $p$ and $n$. Computing $X'X$ is $O(p^2n)$, inverting it is $O(p^3)$, computing $X(X'X)^{-1}$ is $O(np^2)$, and multiplying again by $X'$ is $O(np^2)$ too. So finally the complexity is $O(np^2)$. It is the same for SVD, and I think it is the same for QR (at least if you need only the $p$ first columns of Q). So, who’s doing a benchmark? :) $\endgroup$
    – Elvis
    Commented Mar 3, 2015 at 9:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.