/* generate random square orthogonal matrix via QR decomposition */
static void
test_random_matrix_orth(gsl_matrix *m, const gsl_rng *r)
{
const size_t M = m->size1;
gsl_matrix *A = gsl_matrix_alloc(M, M);
gsl_vector *tau = gsl_vector_alloc(M);
gsl_matrix *R = gsl_matrix_alloc(M, M);
test_random_matrix(A, r, -1.0, 1.0);
gsl_linalg_QR_decomp(A, tau);
gsl_linalg_QR_unpack(A, tau, m, R);
gsl_matrix_free(A);
gsl_matrix_free(R);
gsl_vector_free(tau);
}
/* construct ill-conditioned matrix via SVD */
static void
test_random_matrix_ill(gsl_matrix *m, const gsl_rng *r)
{
const size_t M = m->size1;
const size_t N = m->size2;
gsl_matrix *U = gsl_matrix_alloc(M, M);
gsl_matrix *V = gsl_matrix_alloc(N, N);
gsl_vector *S = gsl_vector_alloc(N);
gsl_matrix_view Uv = gsl_matrix_submatrix(U, 0, 0, M, N);
const double smin = 16.0 * GSL_DBL_EPSILON;
const double smax = 10.0;
const double ratio = pow(smin / smax, 1.0 / (N - 1.0));
double s;
size_t j;
test_random_matrix_orth(U, r);
test_random_matrix_orth(V, r);
/* compute U * S */
s = smax;
for (j = 0; j < N; ++j)
{
gsl_vector_view uj = gsl_matrix_column(U, j);
gsl_vector_scale(&uj.vector, s);
s *= ratio;
}
/* compute m = (U * S) * V' */
gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, &Uv.matrix, V, 0.0, m);
gsl_matrix_free(U);
gsl_matrix_free(V);
gsl_vector_free(S);
}
/* solve system with lambda = 0 and test against OLS solution */
static void
test_reg1(const gsl_matrix * X, const gsl_vector * y,
const gsl_vector * wts, const double tol,
gsl_multifit_linear_workspace * w, const char * desc)
{
const size_t n = X->size1;
const size_t p = X->size2;
double rnorm, snorm, chisq;
gsl_vector *c0 = gsl_vector_alloc(p);
gsl_vector *c1 = gsl_vector_alloc(p);
gsl_matrix *cov = gsl_matrix_alloc(p, p);
size_t j;
if (wts)
{
gsl_matrix *Xs = gsl_matrix_alloc(n, p);
gsl_vector *ys = gsl_vector_alloc(n);
gsl_multifit_wlinear(X, wts, y, c0, cov, &chisq, w);
gsl_multifit_linear_wstdform1(NULL, X, wts, y, Xs, ys, w);
gsl_multifit_linear_svd(Xs, w);
gsl_multifit_linear_solve(0.0, Xs, ys, c1, &rnorm, &snorm, w);
gsl_matrix_free(Xs);
gsl_vector_free(ys);
}
else
{
gsl_multifit_linear(X, y, c0, cov, &chisq, w);
gsl_multifit_linear_svd(X, w);
gsl_multifit_linear_solve(0.0, X, y, c1, &rnorm, &snorm, w);
}
gsl_test_rel(rnorm*rnorm, chisq, tol,
"test_reg1: %s, lambda = 0, n=%zu p=%zu chisq", desc, n, p);
/* test c0 = c1 */
for (j = 0; j < p; ++j)
{
double c0j = gsl_vector_get(c0, j);
double c1j = gsl_vector_get(c1, j);
gsl_test_rel(c1j, c0j, tol, "test_reg1: %s, lambda = 0, n=%zu p=%zu c0/c1",
desc, n, p);
}
gsl_vector_free(c0);
gsl_vector_free(c1);
gsl_matrix_free(cov);
}
/* solve standard form system with given lambda and test against
* normal equations solution, L = I */
static void
test_reg2(const double lambda, const gsl_matrix * X, const gsl_vector * y,
const gsl_vector * wts, const double tol,
gsl_multifit_linear_workspace * w, const char * desc)
{
const size_t n = X->size1;
const size_t p = X->size2;
double rnorm0, snorm0;
double rnorm1, snorm1;
gsl_vector *c0 = gsl_vector_alloc(p);
gsl_vector *c1 = gsl_vector_alloc(p);
gsl_matrix *XTX = gsl_matrix_alloc(p, p); /* X^T W X + lambda^2 I */
gsl_vector *XTy = gsl_vector_alloc(p); /* X^T W y */
gsl_matrix *Xs = gsl_matrix_alloc(n, p);
gsl_vector *ys = gsl_vector_alloc(n);
gsl_vector_view xtx_diag = gsl_matrix_diagonal(XTX);
gsl_permutation *perm = gsl_permutation_alloc(p);
gsl_vector *r = gsl_vector_alloc(n);
int signum;
size_t j;
/* compute Xs = sqrt(W) X and ys = sqrt(W) y */
gsl_multifit_linear_wstdform1(NULL, X, wts, y, Xs, ys, w);
/* construct XTy = X^T W y */
gsl_blas_dgemv(CblasTrans, 1.0, Xs, ys, 0.0, XTy);
/* construct XTX = X^T W X + lambda^2 I */
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, Xs, Xs, 0.0, XTX);
gsl_vector_add_constant(&xtx_diag.vector, lambda*lambda);
/* solve XTX c = XTy with LU decomp */
gsl_linalg_LU_decomp(XTX, perm, &signum);
gsl_linalg_LU_solve(XTX, perm, XTy, c0);
/* compute SVD of X */
gsl_multifit_linear_svd(Xs, w);
/* solve regularized standard form system with lambda */
gsl_multifit_linear_solve(lambda, Xs, ys, c1, &rnorm0, &snorm0, w);
/* test snorm = ||c1|| */
snorm1 = gsl_blas_dnrm2(c1);
gsl_test_rel(snorm0, snorm1, tol, "test_reg2: %s, snorm lambda=%g n=%zu p=%zu",
desc, lambda, n, p);
/* test rnorm = ||y - X c1|| */
gsl_vector_memcpy(r, ys);
gsl_blas_dgemv(CblasNoTrans, -1.0, Xs, c1, 1.0, r);
rnorm1 = gsl_blas_dnrm2(r);
gsl_test_rel(rnorm0, rnorm1, tol, "test_reg2: %s, rnorm lambda=%g n=%zu p=%zu",
desc, lambda, n, p);
/* test c0 = c1 */
for (j = 0; j < p; ++j)
{
double c0j = gsl_vector_get(c0, j);
double c1j = gsl_vector_get(c1, j);
gsl_test_rel(c1j, c0j, tol, "test_reg2: %s, c0/c1 lambda=%g n=%zu p=%zu",
desc, lambda, n, p);
}
gsl_matrix_free(XTX);
gsl_vector_free(XTy);
gsl_matrix_free(Xs);
gsl_vector_free(ys);
gsl_vector_free(c0);
gsl_vector_free(c1);
gsl_vector_free(r);
gsl_permutation_free(perm);
}
/* solve system with given lambda and L = diag(L) and test against
* normal equations solution */
static void
test_reg3(const double lambda, const gsl_vector * L, const gsl_matrix * X,
const gsl_vector * y, const gsl_vector * wts, const double tol,
gsl_multifit_linear_workspace * w, const char * desc)
{
const size_t n = X->size1;
const size_t p = X->size2;
double rnorm0, snorm0;
double rnorm1, snorm1;
gsl_vector *c0 = gsl_vector_alloc(p);
gsl_vector *c1 = gsl_vector_alloc(p);
gsl_matrix *XTX = gsl_matrix_alloc(p, p); /* X^T W X + lambda^2 L^T L */
gsl_vector *XTy = gsl_vector_alloc(p); /* X^T W y */
gsl_matrix *Xs = gsl_matrix_alloc(n, p); /* standard form X~ */
gsl_vector *ys = gsl_vector_alloc(n); /* standard form y~ */
gsl_vector *Lc = gsl_vector_alloc(p);
gsl_vector *r = gsl_vector_alloc(n);
gsl_permutation *perm = gsl_permutation_alloc(p);
int signum;
size_t j;
/* compute Xs = sqrt(W) X, ys = sqrt(W) y */
gsl_multifit_linear_wstdform1(NULL, X, wts, y, Xs, ys, w);
/* construct XTy = X^T W y */
gsl_blas_dgemv(CblasTrans, 1.0, Xs, ys, 0.0, XTy);
/* construct XTX = X^T W X + lambda^2 L^T L */
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, Xs, Xs, 0.0, XTX);
for (j = 0; j < p; ++j)
{
double lj = gsl_vector_get(L, j);
*gsl_matrix_ptr(XTX, j, j) += pow(lambda * lj, 2.0);
}
/* solve XTX c = XTy with LU decomp */
gsl_linalg_LU_decomp(XTX, perm, &signum);
gsl_linalg_LU_solve(XTX, perm, XTy, c0);
/* solve with reg routine */
gsl_multifit_linear_wstdform1(L, X, wts, y, Xs, ys, w);
gsl_multifit_linear_svd(Xs, w);
gsl_multifit_linear_solve(lambda, Xs, ys, c1, &rnorm0, &snorm0, w);
gsl_multifit_linear_genform1(L, c1, c1, w);
/* test snorm = ||L c1|| */
gsl_vector_memcpy(Lc, c1);
gsl_vector_mul(Lc, L);
snorm1 = gsl_blas_dnrm2(Lc);
gsl_test_rel(snorm0, snorm1, tol, "test_reg3: %s, snorm lambda=%g n=%zu p=%zu",
desc, lambda, n, p);
/* test rnorm = ||y - X c1||, compute again Xs = sqrt(W) X and ys = sqrt(W) y */
gsl_multifit_linear_wstdform1(NULL, X, wts, y, Xs, ys, w);
gsl_vector_memcpy(r, ys);
gsl_blas_dgemv(CblasNoTrans, -1.0, Xs, c1, 1.0, r);
rnorm1 = gsl_blas_dnrm2(r);
gsl_test_rel(rnorm0, rnorm1, tol, "test_reg3: %s, rnorm lambda=%g n=%zu p=%zu",
desc, lambda, n, p);
/* test c0 = c1 */
for (j = 0; j < p; ++j)
{
double c0j = gsl_vector_get(c0, j);
double c1j = gsl_vector_get(c1, j);
gsl_test_rel(c1j, c0j, tol, "test_reg3: %s, c0/c1 j=%zu lambda=%g n=%zu p=%zu",
desc, j, lambda, n, p);
}
gsl_matrix_free(Xs);
gsl_matrix_free(XTX);
gsl_vector_free(XTy);
gsl_vector_free(c0);
gsl_vector_free(c1);
gsl_vector_free(Lc);
gsl_vector_free(ys);
gsl_vector_free(r);
gsl_permutation_free(perm);
}
/* solve system with given lambda and L and test against
* normal equations solution */
static void
test_reg4(const double lambda, const gsl_matrix * L, const gsl_matrix * X,
const gsl_vector * y, const gsl_vector * wts, const double tol,
gsl_multifit_linear_workspace * w, const char *desc)
{
const size_t m = L->size1;
const size_t n = X->size1;
const size_t p = X->size2;
double rnorm0, snorm0;
double rnorm1, snorm1;
gsl_vector *c0 = gsl_vector_alloc(p);
gsl_vector *c1 = gsl_vector_alloc(p);
gsl_matrix *LTL = gsl_matrix_alloc(p, p); /* L^T L */
gsl_matrix *XTX = gsl_matrix_alloc(p, p); /* X^T W X + lambda^2 L^T L */
gsl_vector *XTy = gsl_vector_alloc(p); /* X^T W y */
gsl_permutation *perm = gsl_permutation_alloc(p);
gsl_matrix *Xs = (m < p) ? gsl_matrix_alloc(n - (p - m), m) : gsl_matrix_alloc(n, p);
gsl_vector *ys = (m < p) ? gsl_vector_alloc(n - (p - m)) : gsl_vector_alloc(n);
gsl_matrix *M = (m < p) ? gsl_matrix_alloc(n, p) : gsl_matrix_alloc(m, p);
gsl_vector *cs = (m < p) ? gsl_vector_alloc(m) : gsl_vector_alloc(p);
gsl_matrix *WX = gsl_matrix_alloc(n, p);
gsl_vector *Wy = gsl_vector_alloc(n);
gsl_vector *Lc = gsl_vector_alloc(m);
gsl_vector *r = gsl_vector_alloc(n);
gsl_matrix *LQR = gsl_matrix_alloc(m, p);
gsl_vector *Ltau = gsl_vector_alloc(GSL_MIN(m, p));
int signum;
size_t j;
/* compute WX = sqrt(W) X, Wy = sqrt(W) y */
gsl_multifit_linear_wstdform1(NULL, X, wts, y, WX, Wy, w);
/* construct XTy = X^T W y */
gsl_blas_dgemv(CblasTrans, 1.0, WX, Wy, 0.0, XTy);
/* construct XTX = X^T W X + lambda^2 L^T L */
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, L, L, 0.0, LTL);
gsl_matrix_scale(LTL, lambda * lambda);
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, WX, WX, 0.0, XTX);
gsl_matrix_add(XTX, LTL);
/* solve XTX c = XTy with LU decomp */
gsl_linalg_LU_decomp(XTX, perm, &signum);
gsl_linalg_LU_solve(XTX, perm, XTy, c0);
/* solve with reg routine */
gsl_matrix_memcpy(LQR, L);
gsl_multifit_linear_L_decomp(LQR, Ltau);
gsl_multifit_linear_wstdform2(LQR, Ltau, X, wts, y, Xs, ys, M, w);
gsl_multifit_linear_svd(Xs, w);
gsl_multifit_linear_solve(lambda, Xs, ys, cs, &rnorm0, &snorm0, w);
gsl_multifit_linear_wgenform2(LQR, Ltau, X, wts, y, cs, M, c1, w);
/* test snorm = ||L c1|| */
gsl_blas_dgemv(CblasNoTrans, 1.0, L, c1, 0.0, Lc);
snorm1 = gsl_blas_dnrm2(Lc);
gsl_test_rel(snorm0, snorm1, tol, "test_reg4: %s snorm lambda=%g", desc, lambda);
/* test rnorm = ||y - X c1||_W */
gsl_vector_memcpy(r, Wy);
gsl_blas_dgemv(CblasNoTrans, -1.0, WX, c1, 1.0, r);
rnorm1 = gsl_blas_dnrm2(r);
gsl_test_rel(rnorm0, rnorm1, tol, "test_reg4: %s rnorm lambda=%g", desc, lambda);
/* test c0 = c1 */
for (j = 0; j < p; ++j)
{
double c0j = gsl_vector_get(c0, j);
double c1j = gsl_vector_get(c1, j);
gsl_test_rel(c1j, c0j, tol, "test_reg4: %s lambda=%g n=%zu p=%zu j=%zu",
desc, lambda, n, p, j);
}
gsl_matrix_free(LTL);
gsl_matrix_free(XTX);
gsl_vector_free(XTy);
gsl_vector_free(c0);
gsl_vector_free(c1);
gsl_permutation_free(perm);
gsl_matrix_free(Xs);
gsl_vector_free(ys);
gsl_vector_free(cs);
gsl_matrix_free(M);
gsl_vector_free(Lc);
gsl_matrix_free(WX);
gsl_vector_free(Wy);
gsl_vector_free(r);
gsl_matrix_free(LQR);
gsl_vector_free(Ltau);
}
static void
test_reg_system(const size_t n, const size_t p, const gsl_rng *r)
{
gsl_matrix *X = gsl_matrix_alloc(n, p);
gsl_vector *y = gsl_vector_alloc(n);
gsl_vector *c = gsl_vector_alloc(p);
gsl_vector *wts = gsl_vector_alloc(n);
gsl_multifit_linear_workspace *w = gsl_multifit_linear_alloc(n, p);
gsl_multifit_linear_workspace *wbig = gsl_multifit_linear_alloc(n + 10, p + 5);
gsl_vector *diagL = gsl_vector_alloc(p);
gsl_matrix *Lsqr = gsl_matrix_alloc(p, p);
gsl_matrix *Ltall = gsl_matrix_alloc(5*p, p);
gsl_matrix *L1 = gsl_matrix_alloc(p - 1, p);
gsl_matrix *L2 = gsl_matrix_alloc(p - 2, p);
gsl_matrix *L3 = gsl_matrix_alloc(p - 3, p);
gsl_matrix *L5 = gsl_matrix_alloc(p - 5, p);
size_t i;
/* generate random weights */
test_random_vector(wts, r, 0.0, 1.0);
/* generate well-conditioned system and test against OLS solution */
test_random_matrix(X, r, -1.0, 1.0);
test_random_vector(y, r, -1.0, 1.0);
test_reg1(X, y, NULL, 1.0e-10, w, "unweighted");
test_reg1(X, y, wts, 1.0e-10, w, "weighted");
/* generate ill-conditioned system */
test_random_matrix_ill(X, r);
test_random_vector(c, r, -1.0, 1.0);
/* compute y = X c + noise */
gsl_blas_dgemv(CblasNoTrans, 1.0, X, c, 0.0, y);
test_random_vector_noise(r, y);
/* random diag(L) vector */
test_random_vector(diagL, r, -2.0, 2.0);
/* random square and tall L matrices */
test_random_matrix(Lsqr, r, -2.0, 2.0);
test_random_matrix(Ltall, r, -2.0, 2.0);
gsl_multifit_linear_Lk(p, 1, L1);
gsl_multifit_linear_Lk(p, 2, L2);
gsl_multifit_linear_Lk(p, 3, L3);
gsl_multifit_linear_Lk(p, 5, L5);
for (i = 0; i < 3; ++i)
{
/*
* can't make lambda too small or normal equations
* approach won't work well
*/
double lambda = pow(10.0, -(double) i);
/* test unweighted */
test_reg2(lambda, X, y, NULL, 1.0e-6, w, "unweighted");
test_reg3(lambda, diagL, X, y, NULL, 1.0e-6, w, "unweighted");
test_reg4(lambda, Lsqr, X, y, NULL, 1.0e-8, w, "Lsqr unweighted");
test_reg4(lambda, Ltall, X, y, NULL, 1.0e-8, w, "Ltall unweighted");
test_reg4(lambda, L1, X, y, NULL, 1.0e-6, w, "L1 unweighted");
test_reg4(lambda, L2, X, y, NULL, 1.0e-6, w, "L2 unweighted");
test_reg4(lambda, L3, X, y, NULL, 1.0e-5, w, "L3 unweighted");
test_reg4(lambda, L5, X, y, NULL, 1.0e-4, w, "L5 unweighted");
/* test weighted */
test_reg2(lambda, X, y, wts, 1.0e-6, w, "weighted");
test_reg3(lambda, diagL, X, y, wts, 1.0e-6, w, "weighted");
test_reg4(lambda, Lsqr, X, y, wts, 1.0e-8, w, "Lsqr weighted");
test_reg4(lambda, L1, X, y, wts, 1.0e-6, w, "L1 weighted");
test_reg4(lambda, L2, X, y, wts, 1.0e-6, w, "L2 weighted");
test_reg4(lambda, L3, X, y, wts, 1.0e-5, w, "L3 weighted");
test_reg4(lambda, L5, X, y, wts, 1.0e-4, w, "L5 weighted");
/* test again with larger workspace */
test_reg2(lambda, X, y, NULL, 1.0e-6, wbig, "unweighted big");
test_reg3(lambda, diagL, X, y, NULL, 1.0e-6, wbig, "unweighted big");
test_reg4(lambda, Lsqr, X, y, NULL, 1.0e-8, wbig, "Lsqr unweighted big");
test_reg4(lambda, L1, X, y, NULL, 1.0e-6, wbig, "L1 unweighted big");
test_reg4(lambda, L2, X, y, NULL, 1.0e-6, wbig, "L2 unweighted big");
test_reg4(lambda, L3, X, y, NULL, 1.0e-5, wbig, "L3 unweighted big");
test_reg4(lambda, L5, X, y, NULL, 1.0e-4, wbig, "L5 unweighted big");
test_reg2(lambda, X, y, wts, 1.0e-6, wbig, "weighted big");
test_reg3(lambda, diagL, X, y, wts, 1.0e-6, wbig, "weighted big");
test_reg4(lambda, Lsqr, X, y, wts, 1.0e-8, wbig, "Lsqr weighted big");
test_reg4(lambda, L1, X, y, wts, 1.0e-6, wbig, "L1 weighted big");
test_reg4(lambda, L2, X, y, wts, 1.0e-6, wbig, "L2 weighted big");
test_reg4(lambda, L3, X, y, wts, 1.0e-5, wbig, "L3 weighted big");
test_reg4(lambda, L5, X, y, wts, 1.0e-4, wbig, "L5 weighted big");
}
gsl_matrix_free(X);
gsl_vector_free(y);
gsl_vector_free(c);
gsl_vector_free(wts);
gsl_vector_free(diagL);
gsl_matrix_free(Lsqr);
gsl_matrix_free(Ltall);
gsl_matrix_free(L1);
gsl_matrix_free(L2);
gsl_matrix_free(L3);
gsl_matrix_free(L5);
gsl_multifit_linear_free(w);
gsl_multifit_linear_free(wbig);
}
static void
test_reg_sobolev(const size_t p, const size_t kmax, const gsl_rng *r)
{
const double tol = 1.0e-12;
size_t i, j, k;
gsl_matrix *L = gsl_matrix_calloc(p, p);
gsl_matrix *LTL = gsl_matrix_alloc(p, p); /* Sobolov L^T L */
gsl_matrix *LTL2 = gsl_matrix_alloc(p, p); /* alternate L^T L */
gsl_matrix *Li = gsl_matrix_alloc(p, p);
gsl_multifit_linear_workspace *w = gsl_multifit_linear_alloc(p, p);
for (k = 0; k <= kmax; ++k)
{
gsl_vector *alpha = gsl_vector_alloc(k + 1);
/* random weights */
test_random_vector(alpha, r, 0.0, 1.0);
/* compute Sobolev matrix */
gsl_multifit_linear_Lsobolev(p, k, alpha, L, w);
/* compute LTL = L^T L */
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, L, L, 0.0, LTL);
/* now compute LTL2 = L^T L using individual L_i factors */
{
gsl_matrix_set_zero(LTL2);
for (i = 0; i <= k; ++i)
{
gsl_matrix_view Liv = gsl_matrix_submatrix(Li, 0, 0, p - i, p);
double ai = gsl_vector_get(alpha, i);
/* compute a_i L_i */
gsl_multifit_linear_Lk(p, i, &Liv.matrix);
gsl_matrix_scale(&Liv.matrix, ai);
/* LTL += L_i^T L_i */
gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &Liv.matrix, &Liv.matrix, 1.0, LTL2);
}
}
/* test LTL = LTL2 */
for (i = 0; i < p; ++i)
{
for (j = 0; j < p; ++j)
{
double aij = gsl_matrix_get(LTL, i, j);
double bij = gsl_matrix_get(LTL2, i, j);
gsl_test_rel(aij, bij, tol, "sobolov k=%zu LTL(%zu,%zu)", k, i, j);
}
}
gsl_vector_free(alpha);
}
gsl_matrix_free(L);
gsl_matrix_free(Li);
gsl_matrix_free(LTL);
gsl_matrix_free(LTL2);
gsl_multifit_linear_free(w);
}
/* test linear regularized regression */
static void
test_reg(void)
{
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
test_reg_system(100, 15, r);
test_reg_system(100, 50, r);
test_reg_system(100, 99, r);
test_reg_sobolev(20, 10, r);
gsl_rng_free(r);
}