0

I want to use Generalized Cross Validation to compute the optimal ridge parameter for a ridge regression.

The formula due to Golub et al (1979) is

GCV(k)=(1/n Sum_i=1^n e_i(k)^2)/(1/n Sum_i=1^n(1-h_ii(k)))^2

where e_1(k)....e_n(k) are the residuals from applying the ridge regression for a given k.

H(K)=X(X'X+kI)^{-1}X'.

where X is the design matrix.

and h_ii(k) are the digonal elelments of H(k).

I have the ridge regression formula so I can easily compute the numerator.

However I am not sure how to compute the denominator?

I found this set of tools which features a function gcv() but I can't understand how it it calculating the G function?

Indeed it does not appear to be using x in it's calculations? Is this function going what I want above? If so can someone please explain to me how it works?

Bazman
  • 2,058
  • 9
  • 45
  • 65

2 Answers2

2

For certain problems it is possible to compute the generalized cross validation without actually doing any cross validation. For plain (non-regularized) ordinary least squares (OLS, for the searchbots), trace(H_ols), assuming X'X is actually invertible, is = trace(Xinv(X'X)X') = trace(X'X*inv(X'X)) = trace(I) = p which is the number of parameters in the OLS model. This is a constant, for any problem of the same size.

I am not an expert on Ridge Regression, but perhaps gcv.m is doing something similar? After all trace(H_ridge) = trace(X*inv(X.'*X+k*ones(size(X)))*X.') = trace(X'X * inv(X'X + k*ones))

Sorry for the ASCII math.

kousu
  • 79
  • 5
1
(1/n Sum_i=1^n(1-h_ii(k)))^2=(sum(1-diag(H))/n)^2
David
  • 8,449
  • 1
  • 22
  • 32
  • Sorry I should have made it clearer but I am not sure how to calculate H. – Bazman Nov 10 '13 at 22:31
  • Try `H=X*inv(X.'*X+k*ones(size(X)))*X.'` if `k` is treated like a constant. – David Nov 10 '13 at 22:34
  • Hi David, I believe that the inv() in Matlab is not good, see here: http://stackoverflow.com/questions/6509298/is-there-a-fast-way-to-invert-a-matrix-in-matlab if possible I think it's better to use a backslash based solution (sorry i should probably have mode this clear too)!! – Bazman Nov 10 '13 at 22:56
  • That's the problem though I don't know how to forumulate this as a \ type solution. In a pure ridge regression you can append a diagonal matrix with sqrt(k) on the diagonal to X=X^* so that X*X*'= (XX+Ik)^{-1} and a corresponding set of zeros to Y. You can then use the ldivide as per usual. But here X* and X would be of different sizes so it seems using ldivide is not straightforward? – Bazman Nov 10 '13 at 23:33
  • If `X` is m by n, then `X'` is n by m and `X'*X` is n by n and `X*X'` is m by m, so they are square. The link you gave shows how to use backslash to get the inverse, just use `invA = A\eye(size(A));`. – David Nov 10 '13 at 23:37
  • No after the augmenation X=(n,m) and X*=(n+m,m). So X*X*'=(n+m, n+m) but X is still just (n,m). Maybe I need to augement X too (with zeroes?) I'm away from my PC right now so I can't check this. – Bazman Nov 11 '13 at 05:34
  • Oh, usually `X'` denotes the matrix transpose. I don't know about augmentation, sorry. – David Nov 11 '13 at 06:01