Ridge Regression Model with Generalized Cross Validation

An introduction to GCV can be found here:

http://sfb649.wiwi.hu-berlin.de/fedc_homepage/xplore/ebooks/html/csa/node123.html

In ridge regression model with generalized cross validation, we can determine the ridge parameter by selecting \$\lambda\$ that minimizes the following:


where matrix \$H(\lambda)\$ is given as
\$X(X^{T}X+\lambda I)X^{T}\$.

The code:

 void CRidgeRegressionGCV::ComputeRegression()  
 {  
      double minGCV = FLT_MAX;  
      for (int t = 0; t < m_alphas.size(); t++)  
      {  
           double alpha = m_alphas[t];  
           int m = m_X.size();  
           int n = m_X[0].size();  
           MatrixXd mat(m, n + 1);  
           VectorXd rhs(m);  
           for (int i = 0; i < m; i++)  
           {  
                for (int j = 0; j < n + 1; j++)  
                {  
                     if (j < n)  
                     {  
                          mat(i, j) = m_X[i][j];  
                          rhs(i) = m_Y[i];  
                     }  
                     else  
                     {  
                          mat(i, j) = 1;  
                     }  
                }  
           }  
           MatrixXd tmp = mat.transpose()*mat + alpha*MatrixXd::Identity(n + 1, n + 1);  
           JacobiSVD<MatrixXd> svd(tmp, ComputeThinU | ComputeThinV);  
           MatrixXd res = svd.solve(mat.transpose()*rhs);  
           vector<double> coef(n);  
           double intercept;  
           coef.resize(n);  
           for (int i = 0; i < n; i++)  
           {  
                coef[i] = res(i, 0);  
           }  
           intercept = res(n, 0);  
           double GCV = 0;  
           for (int i = 0; i < m; i++)  
           {  
                double sum = 0;  
                for (int j = 0; j < n; j++)  
                {  
                     sum += m_X[i][j] * coef[j];  
                }  
                sum += intercept;  
                double diff = m_Y[i] - sum;  
                double diffsquare = diff*diff;  
                GCV += diffsquare;  
           }  
           GCV /= m;  
           MatrixXd tmp2=svd.solve(mat.transpose());  
           MatrixXd H = mat*tmp2;  
           GCV /= (1-H.trace()/m)*(1-H.trace()/m);  
           if (GCV < minGCV)  
           {  
                minGCV = GCV;  
                m_coef = coef;  
                m_intercept = intercept;  
                m_alpha = alpha;  
           }  
      }  
 }