1 /* multifit/multilinear.c
2 *
3 * Copyright (C) 2000, 2007, 2010 Brian Gough
4 * Copyright (C) 2013, 2015 Patrick Alken
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 #include <config.h>
22 #include <gsl/gsl_errno.h>
23 #include <gsl/gsl_multifit.h>
24 #include <gsl/gsl_blas.h>
25 #include <gsl/gsl_vector.h>
26 #include <gsl/gsl_matrix.h>
27 #include <gsl/gsl_linalg.h>
28
29 #include "linear_common.c"
30
31 static int multifit_linear_svd (const gsl_matrix * X,
32 const int balance,
33 gsl_multifit_linear_workspace * work);
34
35 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)36 gsl_multifit_linear (const gsl_matrix * X,
37 const gsl_vector * y,
38 gsl_vector * c,
39 gsl_matrix * cov,
40 double *chisq, gsl_multifit_linear_workspace * work)
41 {
42 size_t rank;
43 int status = gsl_multifit_linear_tsvd(X, y, GSL_DBL_EPSILON, c, cov, chisq, &rank, work);
44
45 return status;
46 }
47
48 /*
49 gsl_multifit_linear_tsvd()
50 Solve linear least squares system with truncated SVD
51
52 Inputs: X - least squares matrix, n-by-p
53 y - right hand side vector, n-by-1
54 tol - tolerance for singular value truncation; if
55 s_j <= tol * s_0
56 then it is discarded from series expansion
57 c - (output) solution vector, p-by-1
58 cov - (output) covariance matrix, p-by-p
59 chisq - (output) cost function chi^2
60 rank - (output) effective rank (number of singular values
61 used in solution)
62 work - workspace
63 */
64
65 int
gsl_multifit_linear_tsvd(const gsl_matrix * X,const gsl_vector * y,const double tol,gsl_vector * c,gsl_matrix * cov,double * chisq,size_t * rank,gsl_multifit_linear_workspace * work)66 gsl_multifit_linear_tsvd (const gsl_matrix * X,
67 const gsl_vector * y,
68 const double tol,
69 gsl_vector * c,
70 gsl_matrix * cov,
71 double * chisq,
72 size_t * rank,
73 gsl_multifit_linear_workspace * work)
74 {
75 const size_t n = X->size1;
76 const size_t p = X->size2;
77
78 if (y->size != n)
79 {
80 GSL_ERROR("number of observations in y does not match matrix",
81 GSL_EBADLEN);
82 }
83 else if (p != c->size)
84 {
85 GSL_ERROR ("number of parameters c does not match matrix",
86 GSL_EBADLEN);
87 }
88 else if (tol <= 0)
89 {
90 GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
91 }
92 else
93 {
94 int status;
95 double rnorm = 0.0, snorm;
96
97 /* compute balanced SVD */
98 status = gsl_multifit_linear_bsvd (X, work);
99 if (status)
100 return status;
101
102 status = multifit_linear_solve (X, y, tol, -1.0, rank,
103 c, &rnorm, &snorm, work);
104
105 *chisq = rnorm * rnorm;
106
107 /* variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */
108 {
109 double r2 = rnorm * rnorm;
110 double s2 = r2 / (double)(n - *rank);
111 size_t i, j;
112 gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
113 gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);
114
115 for (i = 0; i < p; i++)
116 {
117 gsl_vector_view row_i = gsl_matrix_row (&QSI.matrix, i);
118 double d_i = gsl_vector_get (&D.vector, i);
119
120 for (j = i; j < p; j++)
121 {
122 gsl_vector_view row_j = gsl_matrix_row (&QSI.matrix, j);
123 double d_j = gsl_vector_get (&D.vector, j);
124 double s;
125
126 gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);
127
128 gsl_matrix_set (cov, i, j, s * s2 / (d_i * d_j));
129 gsl_matrix_set (cov, j, i, s * s2 / (d_i * d_j));
130 }
131 }
132 }
133
134 return status;
135 }
136 }
137
138 /*
139 gsl_multifit_linear_svd()
140 Perform SVD decomposition of the matrix X and store in work without
141 balancing
142 */
143
144 int
gsl_multifit_linear_svd(const gsl_matrix * X,gsl_multifit_linear_workspace * work)145 gsl_multifit_linear_svd (const gsl_matrix * X,
146 gsl_multifit_linear_workspace * work)
147 {
148 /* do not balance by default */
149 int status = multifit_linear_svd(X, 0, work);
150
151 return status;
152 }
153
154 /*
155 gsl_multifit_linear_bsvd()
156 Perform SVD decomposition of the matrix X and store in work with
157 balancing
158 */
159
160 int
gsl_multifit_linear_bsvd(const gsl_matrix * X,gsl_multifit_linear_workspace * work)161 gsl_multifit_linear_bsvd (const gsl_matrix * X,
162 gsl_multifit_linear_workspace * work)
163 {
164 int status = multifit_linear_svd(X, 1, work);
165
166 return status;
167 }
168
169 size_t
gsl_multifit_linear_rank(const double tol,const gsl_multifit_linear_workspace * work)170 gsl_multifit_linear_rank(const double tol, const gsl_multifit_linear_workspace * work)
171 {
172 double s0 = gsl_vector_get (work->S, 0);
173 size_t rank = 0;
174 size_t j;
175
176 for (j = 0; j < work->p; j++)
177 {
178 double sj = gsl_vector_get (work->S, j);
179
180 if (sj > tol * s0)
181 ++rank;
182 }
183
184 return rank;
185 }
186
187 /* Estimation of values for given x */
188
189 int
gsl_multifit_linear_est(const gsl_vector * x,const gsl_vector * c,const gsl_matrix * cov,double * y,double * y_err)190 gsl_multifit_linear_est (const gsl_vector * x,
191 const gsl_vector * c,
192 const gsl_matrix * cov, double *y, double *y_err)
193 {
194
195 if (x->size != c->size)
196 {
197 GSL_ERROR ("number of parameters c does not match number of observations x",
198 GSL_EBADLEN);
199 }
200 else if (cov->size1 != cov->size2)
201 {
202 GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
203 }
204 else if (c->size != cov->size1)
205 {
206 GSL_ERROR ("number of parameters c does not match size of covariance matrix cov",
207 GSL_EBADLEN);
208 }
209 else
210 {
211 size_t i, j;
212 double var = 0;
213
214 gsl_blas_ddot(x, c, y); /* y = x.c */
215
216 /* var = x' cov x */
217
218 for (i = 0; i < x->size; i++)
219 {
220 const double xi = gsl_vector_get (x, i);
221 var += xi * xi * gsl_matrix_get (cov, i, i);
222
223 for (j = 0; j < i; j++)
224 {
225 const double xj = gsl_vector_get (x, j);
226 var += 2 * xi * xj * gsl_matrix_get (cov, i, j);
227 }
228 }
229
230 *y_err = sqrt (var);
231
232 return GSL_SUCCESS;
233 }
234 }
235
236 /*
237 gsl_multifit_linear_rcond()
238 Return reciprocal condition number of LS matrix;
239 gsl_multifit_linear_svd() must first be called to
240 compute the SVD of X and its reciprocal condition number
241 */
242
243 double
gsl_multifit_linear_rcond(const gsl_multifit_linear_workspace * w)244 gsl_multifit_linear_rcond (const gsl_multifit_linear_workspace * w)
245 {
246 return w->rcond;
247 }
248
249 /*
250 gsl_multifit_linear_residuals()
251 Compute vector of residuals from fit
252
253 Inputs: X - design matrix
254 y - rhs vector
255 c - fit coefficients
256 r - (output) where to store residuals
257 */
258
259 int
gsl_multifit_linear_residuals(const gsl_matrix * X,const gsl_vector * y,const gsl_vector * c,gsl_vector * r)260 gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
261 const gsl_vector *c, gsl_vector *r)
262 {
263 if (X->size1 != y->size)
264 {
265 GSL_ERROR
266 ("number of observations in y does not match rows of matrix X",
267 GSL_EBADLEN);
268 }
269 else if (X->size2 != c->size)
270 {
271 GSL_ERROR ("number of parameters c does not match columns of matrix X",
272 GSL_EBADLEN);
273 }
274 else if (y->size != r->size)
275 {
276 GSL_ERROR ("number of observations in y does not match number of residuals",
277 GSL_EBADLEN);
278 }
279 else
280 {
281 /* r = y - X c */
282 gsl_vector_memcpy(r, y);
283 gsl_blas_dgemv(CblasNoTrans, -1.0, X, c, 1.0, r);
284
285 return GSL_SUCCESS;
286 }
287 } /* gsl_multifit_linear_residuals() */
288
289 /* Perform a SVD decomposition on the least squares matrix X = U S Q^T
290 *
291 * Inputs: X - least squares matrix
292 * balance - 1 to perform column balancing
293 * work - workspace
294 *
295 * Notes:
296 * 1) On output,
297 * work->A contains the matrix U
298 * work->Q contains the matrix Q
299 * work->S contains the vector of singular values
300 * 2) The matrix X may have smaller dimensions than the workspace
301 * in the case of stdform2() - but the dimensions cannot be larger
302 * 3) On output, work->n and work->p are set to the dimensions of X
303 * 4) On output, work->rcond is set to the reciprocal condition number of X
304 */
305
306 static int
multifit_linear_svd(const gsl_matrix * X,const int balance,gsl_multifit_linear_workspace * work)307 multifit_linear_svd (const gsl_matrix * X,
308 const int balance,
309 gsl_multifit_linear_workspace * work)
310 {
311 const size_t n = X->size1;
312 const size_t p = X->size2;
313
314 if (n > work->nmax || p > work->pmax)
315 {
316 GSL_ERROR("observation matrix larger than workspace", GSL_EBADLEN);
317 }
318 else
319 {
320 gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
321 gsl_matrix_view Q = gsl_matrix_submatrix(work->Q, 0, 0, p, p);
322 gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
323 gsl_vector_view S = gsl_vector_subvector(work->S, 0, p);
324 gsl_vector_view xt = gsl_vector_subvector(work->xt, 0, p);
325 gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);
326
327 /* Copy X to workspace, A <= X */
328
329 gsl_matrix_memcpy (&A.matrix, X);
330
331 /* Balance the columns of the matrix A if requested */
332
333 if (balance)
334 {
335 gsl_linalg_balance_columns (&A.matrix, &D.vector);
336 }
337 else
338 {
339 gsl_vector_set_all (&D.vector, 1.0);
340 }
341
342 /* decompose A into U S Q^T */
343 gsl_linalg_SV_decomp_mod (&A.matrix, &QSI.matrix, &Q.matrix,
344 &S.vector, &xt.vector);
345
346 /* compute reciprocal condition number rcond = smin / smax */
347 {
348 double smin, smax;
349 gsl_vector_minmax(&S.vector, &smin, &smax);
350 work->rcond = smin / smax;
351 }
352
353 work->n = n;
354 work->p = p;
355
356 return GSL_SUCCESS;
357 }
358 }
359