1 /* multifit/multilinear.c
2 *
3 * Copyright (C) 2000, 2007 Brian Gough
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 "gsl__config.h"
21 #include "gsl_errno.h"
22 #include "gsl_multifit.h"
23 #include "gsl_blas.h"
24 #include "gsl_linalg.h"
25
26 /* Fit
27 *
28 * y = X c
29 *
30 * where X is an M x N matrix of M observations for N variables.
31 *
32 */
33
34 int
gsl_multifit_linear(const gsl_matrix * X,const gsl_vector * y,gsl_vector * c,gsl_matrix * cov,double * chisq,gsl_multifit_linear_workspace * work)35 gsl_multifit_linear (const gsl_matrix * X,
36 const gsl_vector * y,
37 gsl_vector * c,
38 gsl_matrix * cov,
39 double *chisq, gsl_multifit_linear_workspace * work)
40 {
41 size_t rank;
42 int status = gsl_multifit_linear_svd (X, y, GSL_DBL_EPSILON, &rank, c,
43 cov, chisq, work);
44 return status;
45 }
46
47 /* Handle the general case of the SVD with tolerance and rank */
48
49 int
gsl_multifit_linear_svd(const gsl_matrix * X,const gsl_vector * y,double tol,size_t * rank,gsl_vector * c,gsl_matrix * cov,double * chisq,gsl_multifit_linear_workspace * work)50 gsl_multifit_linear_svd (const gsl_matrix * X,
51 const gsl_vector * y,
52 double tol,
53 size_t * rank,
54 gsl_vector * c,
55 gsl_matrix * cov,
56 double *chisq, gsl_multifit_linear_workspace * work)
57 {
58 if (X->size1 != y->size)
59 {
60 GSL_ERROR
61 ("number of observations in y does not match rows of matrix X",
62 GSL_EBADLEN);
63 }
64 else if (X->size2 != c->size)
65 {
66 GSL_ERROR ("number of parameters c does not match columns of matrix X",
67 GSL_EBADLEN);
68 }
69 else if (cov->size1 != cov->size2)
70 {
71 GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
72 }
73 else if (c->size != cov->size1)
74 {
75 GSL_ERROR
76 ("number of parameters does not match size of covariance matrix",
77 GSL_EBADLEN);
78 }
79 else if (X->size1 != work->n || X->size2 != work->p)
80 {
81 GSL_ERROR
82 ("size of workspace does not match size of observation matrix",
83 GSL_EBADLEN);
84 }
85 else if (tol <= 0)
86 {
87 GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
88 }
89 else
90 {
91 const size_t n = X->size1;
92 const size_t p = X->size2;
93
94 size_t i, j, p_eff;
95
96 gsl_matrix *A = work->A;
97 gsl_matrix *Q = work->Q;
98 gsl_matrix *QSI = work->QSI;
99 gsl_vector *S = work->S;
100 gsl_vector *xt = work->xt;
101 gsl_vector *D = work->D;
102
103 /* Copy X to workspace, A <= X */
104
105 gsl_matrix_memcpy (A, X);
106
107 /* Balance the columns of the matrix A */
108
109 gsl_linalg_balance_columns (A, D);
110
111 /* Decompose A into U S Q^T */
112
113 gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt);
114
115 /* Solve y = A c for c */
116
117 gsl_blas_dgemv (CblasTrans, 1.0, A, y, 0.0, xt);
118
119 /* Scale the matrix Q, Q' = Q S^-1 */
120
121 gsl_matrix_memcpy (QSI, Q);
122
123 {
124 double alpha0 = gsl_vector_get (S, 0);
125 p_eff = 0;
126
127 for (j = 0; j < p; j++)
128 {
129 gsl_vector_view column = gsl_matrix_column (QSI, j);
130 double alpha = gsl_vector_get (S, j);
131
132 if (alpha <= tol * alpha0) {
133 alpha = 0.0;
134 } else {
135 alpha = 1.0 / alpha;
136 p_eff++;
137 }
138
139 gsl_vector_scale (&column.vector, alpha);
140 }
141
142 *rank = p_eff;
143 }
144
145 gsl_vector_set_zero (c);
146
147 gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c);
148
149 /* Unscale the balancing factors */
150
151 gsl_vector_div (c, D);
152
153 /* Compute chisq, from residual r = y - X c */
154
155 {
156 double s2 = 0, r2 = 0;
157
158 for (i = 0; i < n; i++)
159 {
160 double yi = gsl_vector_get (y, i);
161 gsl_vector_const_view row = gsl_matrix_const_row (X, i);
162 double y_est, ri;
163 gsl_blas_ddot (&row.vector, c, &y_est);
164 ri = yi - y_est;
165 r2 += ri * ri;
166 }
167
168 s2 = r2 / (n - p_eff); /* p_eff == rank */
169
170 *chisq = r2;
171
172 /* Form variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */
173
174 for (i = 0; i < p; i++)
175 {
176 gsl_vector_view row_i = gsl_matrix_row (QSI, i);
177 double d_i = gsl_vector_get (D, i);
178
179 for (j = i; j < p; j++)
180 {
181 gsl_vector_view row_j = gsl_matrix_row (QSI, j);
182 double d_j = gsl_vector_get (D, j);
183 double s;
184
185 gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);
186
187 gsl_matrix_set (cov, i, j, s * s2 / (d_i * d_j));
188 gsl_matrix_set (cov, j, i, s * s2 / (d_i * d_j));
189 }
190 }
191 }
192
193 return GSL_SUCCESS;
194 }
195 }
196
197 int
gsl_multifit_wlinear(const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,gsl_vector * c,gsl_matrix * cov,double * chisq,gsl_multifit_linear_workspace * work)198 gsl_multifit_wlinear (const gsl_matrix * X,
199 const gsl_vector * w,
200 const gsl_vector * y,
201 gsl_vector * c,
202 gsl_matrix * cov,
203 double *chisq, gsl_multifit_linear_workspace * work)
204 {
205 size_t rank;
206 int status = gsl_multifit_wlinear_svd (X, w, y, GSL_DBL_EPSILON, &rank, c,
207 cov, chisq, work);
208 return status;
209 }
210
211
212 int
gsl_multifit_wlinear_svd(const gsl_matrix * X,const gsl_vector * w,const gsl_vector * y,double tol,size_t * rank,gsl_vector * c,gsl_matrix * cov,double * chisq,gsl_multifit_linear_workspace * work)213 gsl_multifit_wlinear_svd (const gsl_matrix * X,
214 const gsl_vector * w,
215 const gsl_vector * y,
216 double tol,
217 size_t * rank,
218 gsl_vector * c,
219 gsl_matrix * cov,
220 double *chisq, gsl_multifit_linear_workspace * work)
221 {
222 if (X->size1 != y->size)
223 {
224 GSL_ERROR
225 ("number of observations in y does not match rows of matrix X",
226 GSL_EBADLEN);
227 }
228 else if (X->size2 != c->size)
229 {
230 GSL_ERROR ("number of parameters c does not match columns of matrix X",
231 GSL_EBADLEN);
232 }
233 else if (w->size != y->size)
234 {
235 GSL_ERROR ("number of weights does not match number of observations",
236 GSL_EBADLEN);
237 }
238 else if (cov->size1 != cov->size2)
239 {
240 GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
241 }
242 else if (c->size != cov->size1)
243 {
244 GSL_ERROR
245 ("number of parameters does not match size of covariance matrix",
246 GSL_EBADLEN);
247 }
248 else if (X->size1 != work->n || X->size2 != work->p)
249 {
250 GSL_ERROR
251 ("size of workspace does not match size of observation matrix",
252 GSL_EBADLEN);
253 }
254 else
255 {
256 const size_t n = X->size1;
257 const size_t p = X->size2;
258
259 size_t i, j, p_eff;
260
261 gsl_matrix *A = work->A;
262 gsl_matrix *Q = work->Q;
263 gsl_matrix *QSI = work->QSI;
264 gsl_vector *S = work->S;
265 gsl_vector *t = work->t;
266 gsl_vector *xt = work->xt;
267 gsl_vector *D = work->D;
268
269 /* Scale X, A = sqrt(w) X */
270
271 gsl_matrix_memcpy (A, X);
272
273 for (i = 0; i < n; i++)
274 {
275 double wi = gsl_vector_get (w, i);
276
277 if (wi < 0)
278 wi = 0;
279
280 {
281 gsl_vector_view row = gsl_matrix_row (A, i);
282 gsl_vector_scale (&row.vector, sqrt (wi));
283 }
284 }
285
286 /* Balance the columns of the matrix A */
287
288 gsl_linalg_balance_columns (A, D);
289
290 /* Decompose A into U S Q^T */
291
292 gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt);
293
294 /* Solve sqrt(w) y = A c for c, by first computing t = sqrt(w) y */
295
296 for (i = 0; i < n; i++)
297 {
298 double wi = gsl_vector_get (w, i);
299 double yi = gsl_vector_get (y, i);
300 if (wi < 0)
301 wi = 0;
302 gsl_vector_set (t, i, sqrt (wi) * yi);
303 }
304
305 gsl_blas_dgemv (CblasTrans, 1.0, A, t, 0.0, xt);
306
307 /* Scale the matrix Q, Q' = Q S^-1 */
308
309 gsl_matrix_memcpy (QSI, Q);
310
311 {
312 double alpha0 = gsl_vector_get (S, 0);
313 p_eff = 0;
314
315 for (j = 0; j < p; j++)
316 {
317 gsl_vector_view column = gsl_matrix_column (QSI, j);
318 double alpha = gsl_vector_get (S, j);
319
320 if (alpha <= tol * alpha0) {
321 alpha = 0.0;
322 } else {
323 alpha = 1.0 / alpha;
324 p_eff++;
325 }
326
327 gsl_vector_scale (&column.vector, alpha);
328 }
329
330 *rank = p_eff;
331 }
332
333 gsl_vector_set_zero (c);
334
335 /* Solution */
336
337 gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c);
338
339 /* Unscale the balancing factors */
340
341 gsl_vector_div (c, D);
342
343 /* Form covariance matrix cov = (Q S^-1) (Q S^-1)^T */
344
345 for (i = 0; i < p; i++)
346 {
347 gsl_vector_view row_i = gsl_matrix_row (QSI, i);
348 double d_i = gsl_vector_get (D, i);
349
350 for (j = i; j < p; j++)
351 {
352 gsl_vector_view row_j = gsl_matrix_row (QSI, j);
353 double d_j = gsl_vector_get (D, j);
354 double s;
355
356 gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);
357
358 gsl_matrix_set (cov, i, j, s / (d_i * d_j));
359 gsl_matrix_set (cov, j, i, s / (d_i * d_j));
360 }
361 }
362
363 /* Compute chisq, from residual r = y - X c */
364
365 {
366 double r2 = 0;
367
368 for (i = 0; i < n; i++)
369 {
370 double yi = gsl_vector_get (y, i);
371 double wi = gsl_vector_get (w, i);
372 gsl_vector_const_view row = gsl_matrix_const_row (X, i);
373 double y_est, ri;
374 gsl_blas_ddot (&row.vector, c, &y_est);
375 ri = yi - y_est;
376 r2 += wi * ri * ri;
377 }
378
379 *chisq = r2;
380 }
381
382 return GSL_SUCCESS;
383 }
384 }
385
386
387 int
gsl_multifit_linear_est(const gsl_vector * x,const gsl_vector * c,const gsl_matrix * cov,double * y,double * y_err)388 gsl_multifit_linear_est (const gsl_vector * x,
389 const gsl_vector * c,
390 const gsl_matrix * cov, double *y, double *y_err)
391 {
392
393 if (x->size != c->size)
394 {
395 GSL_ERROR ("number of parameters c does not match number of observations x",
396 GSL_EBADLEN);
397 }
398 else if (cov->size1 != cov->size2)
399 {
400 GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
401 }
402 else if (c->size != cov->size1)
403 {
404 GSL_ERROR ("number of parameters c does not match size of covariance matrix cov",
405 GSL_EBADLEN);
406 }
407 else
408 {
409 size_t i, j;
410 double var = 0;
411
412 gsl_blas_ddot(x, c, y); /* y = x.c */
413
414 /* var = x' cov x */
415
416 for (i = 0; i < x->size; i++)
417 {
418 const double xi = gsl_vector_get (x, i);
419 var += xi * xi * gsl_matrix_get (cov, i, i);
420
421 for (j = 0; j < i; j++)
422 {
423 const double xj = gsl_vector_get (x, j);
424 var += 2 * xi * xj * gsl_matrix_get (cov, i, j);
425 }
426 }
427
428 *y_err = sqrt (var);
429
430 return GSL_SUCCESS;
431 }
432 }
433