- //Set m-rank(r) singular values to zero, and recompose
- //best rank(r) approximation of the initial product
- for (int j=rank[r]; j<m; j++)
- S->data[j] = 0.0;
-
- //[intermediate step] Compute hatBetaR = U * S * t(V)
- double* U = matrixM->data; //GSL require double precision
- for (int j=0; j<p; j++)
- {
- for (int jj=0; jj<m; jj++)
- {
- Real dotProduct = 0.0;
- for (int u=0; u<m; u++)
- dotProduct += U[j*m+u] * S->data[u] * V->data[jj*m+u];
- hatBetaR[mi(j,jj,p,m)] = dotProduct;
- }
- }
+ //U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
+ if (p >= m)
+ gsl_linalg_SV_decomp(matrixM, V, S, work);
+ else
+ {
+ gsl_matrix* matrixM_ = gsl_matrix_alloc(m, p);
+ for (int j=0; j<m; j++)
+ {
+ for (int jj=0; jj<p; jj++)
+ matrixM_->data[jj*m+j] = matrixM->data[j*p +jj];
+ }
+ gsl_matrix* V_ = gsl_matrix_alloc(p,p);
+ gsl_vector* S_ = gsl_vector_alloc(p);
+ gsl_vector* work_ = gsl_vector_alloc(p);
+ gsl_linalg_SV_decomp(matrixM_, V_, S_, work_);
+ gsl_vector_free(work_);
+ for (int j=0; j<m; j++)
+ {
+ if (j < p)
+ S->data[j] = S_->data[j];
+ else
+ S->data[j] = 0.0;
+ for (int jj=0; jj<m; jj++)
+ {
+ if (j < p && jj < p)
+ V->data[jj * m + j] = V_->data[jj * m + j];
+ else
+ V->data[jj * m + j] = 0.0;
+ }
+ }
+ for (int j=0; j<m; j++)
+ {
+ for (int jj=0; jj<p; jj++)
+ matrixM->data[j*p+jj] = matrixM_->data[jj*m +j];
+ }
+ gsl_matrix_free(matrixM_);
+ gsl_matrix_free(V_);
+ gsl_vector_free(S_);
+ }