1 /* multilarge.c
2 *
3 * Copyright (C) 2015 Patrick Alken
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <gsl/gsl_math.h>
22 #include <gsl/gsl_vector.h>
23 #include <gsl/gsl_matrix.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_multifit.h>
26 #include <gsl/gsl_multilarge.h>
27 #include <gsl/gsl_blas.h>
28
29 gsl_multilarge_linear_workspace *
gsl_multilarge_linear_alloc(const gsl_multilarge_linear_type * T,const size_t p)30 gsl_multilarge_linear_alloc(const gsl_multilarge_linear_type *T,
31 const size_t p)
32 {
33 gsl_multilarge_linear_workspace *w;
34
35 w = calloc(1, sizeof(gsl_multilarge_linear_workspace));
36 if (w == NULL)
37 {
38 GSL_ERROR_NULL("failed to allocate space for workspace",
39 GSL_ENOMEM);
40 }
41
42 w->type = T;
43
44 w->state = w->type->alloc(p);
45 if (w->state == NULL)
46 {
47 gsl_multilarge_linear_free(w);
48 GSL_ERROR_NULL("failed to allocate space for multilarge state",
49 GSL_ENOMEM);
50 }
51
52 w->p = p;
53
54 /* initialize newly allocated state */
55 gsl_multilarge_linear_reset(w);
56
57 return w;
58 }
59
60 void
gsl_multilarge_linear_free(gsl_multilarge_linear_workspace * w)61 gsl_multilarge_linear_free(gsl_multilarge_linear_workspace *w)
62 {
63 RETURN_IF_NULL(w);
64
65 if (w->state)
66 w->type->free(w->state);
67
68 free(w);
69 }
70
71 const char *
gsl_multilarge_linear_name(const gsl_multilarge_linear_workspace * w)72 gsl_multilarge_linear_name(const gsl_multilarge_linear_workspace *w)
73 {
74 return w->type->name;
75 }
76
77 int
gsl_multilarge_linear_reset(gsl_multilarge_linear_workspace * w)78 gsl_multilarge_linear_reset(gsl_multilarge_linear_workspace *w)
79 {
80 int status = w->type->reset(w->state);
81 return status;
82 }
83
84 int
gsl_multilarge_linear_accumulate(gsl_matrix * X,gsl_vector * y,gsl_multilarge_linear_workspace * w)85 gsl_multilarge_linear_accumulate(gsl_matrix * X, gsl_vector * y,
86 gsl_multilarge_linear_workspace * w)
87 {
88 int status = w->type->accumulate(X, y, w->state);
89 return status;
90 }
91
92 int
gsl_multilarge_linear_solve(const double lambda,gsl_vector * c,double * rnorm,double * snorm,gsl_multilarge_linear_workspace * w)93 gsl_multilarge_linear_solve(const double lambda, gsl_vector * c,
94 double * rnorm, double * snorm,
95 gsl_multilarge_linear_workspace * w)
96 {
97 int status = w->type->solve(lambda, c, rnorm, snorm, w->state);
98 return status;
99 }
100
101 int
gsl_multilarge_linear_rcond(double * rcond,gsl_multilarge_linear_workspace * w)102 gsl_multilarge_linear_rcond(double *rcond, gsl_multilarge_linear_workspace * w)
103 {
104 int status = w->type->rcond(rcond, w->state);
105 return status;
106 }
107
108 int
gsl_multilarge_linear_lcurve(gsl_vector * reg_param,gsl_vector * rho,gsl_vector * eta,gsl_multilarge_linear_workspace * w)109 gsl_multilarge_linear_lcurve(gsl_vector * reg_param, gsl_vector * rho,
110 gsl_vector * eta,
111 gsl_multilarge_linear_workspace * w)
112 {
113 const size_t len = reg_param->size;
114
115 if (len != rho->size)
116 {
117 GSL_ERROR ("reg_param and rho have different sizes", GSL_EBADLEN);
118 }
119 else if (len != eta->size)
120 {
121 GSL_ERROR ("reg_param and eta have different sizes", GSL_EBADLEN);
122 }
123 else
124 {
125 int status = w->type->lcurve(reg_param, rho, eta, w->state);
126 return status;
127 }
128 }
129
130 /*
131 gsl_multilarge_linear_wstdform1()
132 Using regularization matrix
133 L = diag(l_1,l_2,...,l_p), transform to Tikhonov standard form:
134
135 X~ = sqrt(W) X L^{-1}
136 y~ = sqrt(W) y
137 c~ = L c
138
139 Inputs: L - Tikhonov matrix as a vector of diagonal elements p-by-1;
140 or NULL for L = I
141 X - least squares matrix n-by-p
142 y - right hand side vector n-by-1
143 w - weight vector n-by-1; or NULL for W = I
144 Xs - least squares matrix in standard form X~ n-by-p
145 ys - right hand side vector in standard form y~ n-by-1
146 work - workspace
147
148 Return: success/error
149
150 Notes:
151 1) It is allowed for X = Xs and y = ys
152 */
153
154 int
gsl_multilarge_linear_wstdform1(const gsl_vector * L,const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_multilarge_linear_workspace * work)155 gsl_multilarge_linear_wstdform1 (const gsl_vector * L,
156 const gsl_matrix * X,
157 const gsl_vector * w,
158 const gsl_vector * y,
159 gsl_matrix * Xs,
160 gsl_vector * ys,
161 gsl_multilarge_linear_workspace * work)
162 {
163 const size_t n = X->size1;
164 const size_t p = X->size2;
165
166 (void) work;
167
168 if (L != NULL && p != L->size)
169 {
170 GSL_ERROR("L vector does not match X", GSL_EBADLEN);
171 }
172 else if (n != y->size)
173 {
174 GSL_ERROR("y vector does not match X", GSL_EBADLEN);
175 }
176 else if (w != NULL && n != w->size)
177 {
178 GSL_ERROR("weight vector does not match X", GSL_EBADLEN);
179 }
180 else if (n != Xs->size1 || p != Xs->size2)
181 {
182 GSL_ERROR("Xs matrix dimensions do not match X", GSL_EBADLEN);
183 }
184 else if (n != ys->size)
185 {
186 GSL_ERROR("ys vector must be length n", GSL_EBADLEN);
187 }
188 else
189 {
190 int status = GSL_SUCCESS;
191
192 /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
193 status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
194 if (status)
195 return status;
196
197 if (L != NULL)
198 {
199 size_t j;
200
201 /* construct X~ = sqrt(W) X * L^{-1} matrix */
202 for (j = 0; j < p; ++j)
203 {
204 gsl_vector_view Xj = gsl_matrix_column(Xs, j);
205 double lj = gsl_vector_get(L, j);
206
207 if (lj == 0.0)
208 {
209 GSL_ERROR("L matrix is singular", GSL_EDOM);
210 }
211
212 gsl_vector_scale(&Xj.vector, 1.0 / lj);
213 }
214 }
215
216 return status;
217 }
218 }
219
220 int
gsl_multilarge_linear_stdform1(const gsl_vector * L,const gsl_matrix * X,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_multilarge_linear_workspace * work)221 gsl_multilarge_linear_stdform1 (const gsl_vector * L,
222 const gsl_matrix * X,
223 const gsl_vector * y,
224 gsl_matrix * Xs,
225 gsl_vector * ys,
226 gsl_multilarge_linear_workspace * work)
227 {
228 int status;
229
230 status = gsl_multilarge_linear_wstdform1(L, X, NULL, y, Xs, ys, work);
231
232 return status;
233 }
234
235 int
gsl_multilarge_linear_L_decomp(gsl_matrix * L,gsl_vector * tau)236 gsl_multilarge_linear_L_decomp (gsl_matrix * L, gsl_vector * tau)
237 {
238 const size_t m = L->size1;
239 const size_t p = L->size2;
240
241 if (m < p)
242 {
243 GSL_ERROR("m < p not yet supported", GSL_EBADLEN);
244 }
245 else
246 {
247 int status;
248
249 status = gsl_multifit_linear_L_decomp(L, tau);
250
251 return status;
252 }
253 }
254
255 int
gsl_multilarge_linear_wstdform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_multilarge_linear_workspace * work)256 gsl_multilarge_linear_wstdform2 (const gsl_matrix * LQR,
257 const gsl_vector * Ltau,
258 const gsl_matrix * X,
259 const gsl_vector * w,
260 const gsl_vector * y,
261 gsl_matrix * Xs,
262 gsl_vector * ys,
263 gsl_multilarge_linear_workspace * work)
264 {
265 const size_t m = LQR->size1;
266 const size_t n = X->size1;
267 const size_t p = X->size2;
268
269 (void) Ltau;
270
271 if (p != work->p)
272 {
273 GSL_ERROR("X has wrong number of columns", GSL_EBADLEN);
274 }
275 else if (p != LQR->size2)
276 {
277 GSL_ERROR("LQR and X matrices have different numbers of columns", GSL_EBADLEN);
278 }
279 else if (n != y->size)
280 {
281 GSL_ERROR("y vector does not match X", GSL_EBADLEN);
282 }
283 else if (w != NULL && n != w->size)
284 {
285 GSL_ERROR("weights vector must be length n", GSL_EBADLEN);
286 }
287 else if (m < p)
288 {
289 GSL_ERROR("m < p not yet supported", GSL_EBADLEN);
290 }
291 else if (n != Xs->size1 || p != Xs->size2)
292 {
293 GSL_ERROR("Xs matrix must be n-by-p", GSL_EBADLEN);
294 }
295 else if (n != ys->size)
296 {
297 GSL_ERROR("ys vector must have length n", GSL_EBADLEN);
298 }
299 else
300 {
301 int status;
302 size_t i;
303 gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p);
304
305 /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
306 status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
307 if (status)
308 return status;
309
310 /* compute X~ = X R^{-1} using QR decomposition of L */
311 for (i = 0; i < n; ++i)
312 {
313 gsl_vector_view v = gsl_matrix_row(Xs, i);
314
315 /* solve: R^T y = X_i */
316 gsl_blas_dtrsv(CblasUpper, CblasTrans, CblasNonUnit, &R.matrix, &v.vector);
317 }
318
319 return GSL_SUCCESS;
320 }
321 }
322
323 int
gsl_multilarge_linear_stdform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_matrix * X,const gsl_vector * y,gsl_matrix * Xs,gsl_vector * ys,gsl_multilarge_linear_workspace * work)324 gsl_multilarge_linear_stdform2 (const gsl_matrix * LQR,
325 const gsl_vector * Ltau,
326 const gsl_matrix * X,
327 const gsl_vector * y,
328 gsl_matrix * Xs,
329 gsl_vector * ys,
330 gsl_multilarge_linear_workspace * work)
331 {
332 int status;
333
334 status = gsl_multilarge_linear_wstdform2(LQR, Ltau, X, NULL, y, Xs, ys, work);
335
336 return status;
337 }
338
339 /*
340 gsl_multilarge_linear_genform1()
341 Backtransform regularized solution vector using matrix
342 L = diag(L)
343 */
344
345 int
gsl_multilarge_linear_genform1(const gsl_vector * L,const gsl_vector * cs,gsl_vector * c,gsl_multilarge_linear_workspace * work)346 gsl_multilarge_linear_genform1 (const gsl_vector * L,
347 const gsl_vector * cs,
348 gsl_vector * c,
349 gsl_multilarge_linear_workspace * work)
350 {
351 if (L->size != work->p)
352 {
353 GSL_ERROR("L vector does not match workspace", GSL_EBADLEN);
354 }
355 else if (L->size != cs->size)
356 {
357 GSL_ERROR("cs vector does not match L", GSL_EBADLEN);
358 }
359 else if (L->size != c->size)
360 {
361 GSL_ERROR("c vector does not match L", GSL_EBADLEN);
362 }
363 else
364 {
365 /* compute true solution vector c = L^{-1} c~ */
366 gsl_vector_memcpy(c, cs);
367 gsl_vector_div(c, L);
368
369 return GSL_SUCCESS;
370 }
371 }
372
373 int
gsl_multilarge_linear_genform2(const gsl_matrix * LQR,const gsl_vector * Ltau,const gsl_vector * cs,gsl_vector * c,gsl_multilarge_linear_workspace * work)374 gsl_multilarge_linear_genform2 (const gsl_matrix * LQR,
375 const gsl_vector * Ltau,
376 const gsl_vector * cs,
377 gsl_vector * c,
378 gsl_multilarge_linear_workspace * work)
379 {
380 const size_t m = LQR->size1;
381 const size_t p = LQR->size2;
382
383 (void) Ltau;
384 (void) work;
385
386 if (p != c->size)
387 {
388 GSL_ERROR("c vector does not match LQR", GSL_EBADLEN);
389 }
390 else if (m < p)
391 {
392 GSL_ERROR("m < p not yet supported", GSL_EBADLEN);
393 }
394 else if (p != cs->size)
395 {
396 GSL_ERROR("cs vector size does not match c", GSL_EBADLEN);
397 }
398 else
399 {
400 int s;
401 gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p); /* R factor of L */
402
403 /* solve R c = cs for true solution c, using QR decomposition of L */
404 gsl_vector_memcpy(c, cs);
405 s = gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &R.matrix, c);
406
407 return s;
408 }
409 }
410
411 const gsl_matrix *
gsl_multilarge_linear_matrix_ptr(const gsl_multilarge_linear_workspace * work)412 gsl_multilarge_linear_matrix_ptr (const gsl_multilarge_linear_workspace * work)
413 {
414 return work->type->matrix_ptr(work->state);
415 }
416
417 const gsl_vector *
gsl_multilarge_linear_rhs_ptr(const gsl_multilarge_linear_workspace * work)418 gsl_multilarge_linear_rhs_ptr (const gsl_multilarge_linear_workspace * work)
419 {
420 return work->type->rhs_ptr(work->state);
421 }
422