single-multivariate.c
1 //! @file single-multivariate.c
2 //! @author J. Marcel van der Veer
3
4 //! @section Copyright
5 //!
6 //! This file is part of Algol68G - an Algol 68 compiler-interpreter.
7 //! Copyright 2001-2023 J. Marcel van der Veer [algol68g@xs4all.nl].
8
9 //! @section License
10 //!
11 //! This program is free software; you can redistribute it and/or modify it
12 //! under the terms of the GNU General Public License as published by the
13 //! Free Software Foundation; either version 3 of the License, or
14 //! (at your option) any later version.
15 //!
16 //! This program is distributed in the hope that it will be useful, but
17 //! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 //! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
19 //! more details. You should have received a copy of the GNU General Public
20 //! License along with this program. If not, see [http://www.gnu.org/licenses/].
21
22 //! @section Synopsis
23 //!
24 //! REAL multivariate regression.
25
26 // This code implements least square regression methods:
27 // OLS - Ordinary Least Squares
28 // PCR - Principle Component Regression, OLS after dimension reduction
29 // TLS - Total Least Squares
30 // PLS - Partial Least Squares, TLS after dimension reduction
31
32 #include "a68g.h"
33
34 #if defined (HAVE_GSL)
35
36 #include "a68g-torrix.h"
37 #include "a68g-prelude-gsl.h"
38
39 #define SMALL_EIGEN 1e-9
40
41 //! @brief column-centered data matrix.
42
43 gsl_matrix *column_mean (gsl_matrix *DATA)
44 {
45 // A is MxN; M samples x N features.
46 unt M = SIZE1 (DATA), N = SIZE2 (DATA);
47 gsl_matrix *C = gsl_matrix_calloc (M, N);
48 for (int i = 0; i < N; i++) {
49 DOUBLE_T sum = 0;
50 for (int j = 0; j < M; j++) {
51 sum += gsl_matrix_get (DATA, j, i);
52 }
53 REAL_T mean = sum / M;
54 for (int j = 0; j < M; j++) {
55 gsl_matrix_set (C, j, i, gsl_matrix_get (DATA, j, i) - mean);
56 }
57 }
58 return C;
59 }
60
61 //! @brief left columns from matrix.
62
63 static gsl_matrix *left_columns (NODE_T *p, gsl_matrix *X, int cols)
64 {
65 MATH_RTE (p, cols < 0, M_INT, "invalid number of columns");
66 unt M = SIZE1 (X), N = SIZE2 (X);
67 if (cols == 0 || cols > N) {
68 cols = N;
69 }
70 gsl_matrix *Y = gsl_matrix_calloc (M, cols);
71 for (int i = 0; i < M; i++) {
72 for (int j = 0; j < cols; j++) {
73 gsl_matrix_set (Y, i, j, gsl_matrix_get (X, i, j));
74 }
75 }
76 return Y;
77 }
78
79 //! @brief Compute Moore-Penrose pseudo inverse.
80
81 void compute_pseudo_inverse (NODE_T *p, gsl_matrix **mpinv, gsl_matrix *X, REAL_T lim)
82 {
83 // The Moore-Penrose pseudo inverse gives a least-square approximate solution
84 // for a system of linear equations not having an exact solution.
85 // Multivariate statistics is a well known application.
86 MATH_RTE (p, X == NO_REAL_MATRIX, M_ROW_ROW_REAL, "empty data matrix");
87 MATH_RTE (p, lim < 0, M_REAL, "invalid limit");
88 unt M = SIZE1 (X), N = SIZE2 (X);
89 // GSL only handles M >= N, transpose commutes with pseudo inverse.
90 gsl_matrix *U;
91 BOOL_T transpose = (M < N);
92 if (transpose) {
93 M = SIZE2 (X); N = SIZE1 (X);
94 U = gsl_matrix_calloc (M, N);
95 gsl_matrix_transpose_memcpy (U, X);
96 } else {
97 U = gsl_matrix_calloc (M, N);
98 gsl_matrix_memcpy (U, X);
99 }
100 // A = USVᵀ by Jacobi, more precise than Golub-Reinsch.
101 // The GSL routine yields V, not Vᵀ. U is decomposed in place.
102 gsl_matrix *V = gsl_matrix_calloc (N, N);
103 gsl_vector *Sv = gsl_vector_calloc (N);
104 ASSERT_GSL (gsl_linalg_SV_decomp_jacobi (U, V, Sv));
105 // Compute S⁻¹.
106 // Very small eigenvalues wreak havoc on a pseudo inverse.
107 // So diagonal elements < 'lim * largest element' are set to zero.
108 // Python NumPy default is 1e-15, but assumes a Hermitian matrix.
109 // SVD gives singular values sorted in descending order.
110 REAL_T x0 = gsl_vector_get (Sv, 0);
111 MATH_RTE (p, x0 == 0, M_ROW_ROW_REAL, "zero eigenvalue or singular value");
112 gsl_matrix *S_inv = gsl_matrix_calloc (N, M);
113 gsl_matrix_set (S_inv, 0, 0, 1 / x0);
114 for (int i = 1; i < N; i++) {
115 REAL_T x = gsl_vector_get (Sv, i);
116 if ((x / x0) > lim) {
117 gsl_matrix_set (S_inv, i, i, 1 / x);
118 } else {
119 gsl_matrix_set (S_inv, i, i, 0);
120 }
121 }
122 a68_vector_free (Sv);
123 // GSL SVD yields thin SVD - pad with zeros.
124 gsl_matrix *Uf = gsl_matrix_calloc (M, M);
125 for (int i = 0; i < M; i++) {
126 for (int j = 0; j < N; j++) {
127 gsl_matrix_set (Uf, i, j, gsl_matrix_get (U, i, j));
128 }
129 }
130 // Compute pseudo inverse A⁻¹ = VS⁻¹Uᵀ.
131 gsl_matrix *VS_inv = NO_REAL_MATRIX, *X_inv = NO_REAL_MATRIX;
132 a68_dgemm (SELF, SELF, 1, V, S_inv, 0, &VS_inv);
133 a68_dgemm (SELF, FLIP, 1, VS_inv, Uf, 0, &X_inv);
134 // Compose result.
135 if (transpose) {
136 (*mpinv) = gsl_matrix_calloc (M, N);
137 gsl_matrix_transpose_memcpy ((*mpinv), X_inv);
138 } else {
139 (*mpinv) = gsl_matrix_calloc (N, M);
140 gsl_matrix_memcpy ((*mpinv), X_inv);
141 }
142 // Clean up.
143 a68_matrix_free (S_inv);
144 a68_matrix_free (U);
145 a68_matrix_free (Uf);
146 a68_matrix_free (V);
147 a68_matrix_free (VS_inv);
148 a68_matrix_free (X_inv);
149 }
150
151 //! @brief Compute Principal Component Analysis by SVD.
152
153 gsl_matrix *pca_svd (gsl_matrix *X, gsl_vector ** Sv)
154 {
155 // In PCA, SVD is closely related to eigenvector decomposition of the covariance matrix:
156 // If Cov = XᵀX = VLVᵀ and X = USVᵀ then XᵀX = VSUᵀUSVᵀ = VS²Vᵀ, hence L = S².
157 // M samples, N features.
158 unt M = SIZE1 (X), N = SIZE2 (X);
159 // GSL does thin SVD, only handles M >= N, overdetermined systems.
160 BOOL_T transpose = (M < N);
161 gsl_matrix *U = NO_REAL_MATRIX;
162 if (transpose) {
163 // Xᵀ = VSUᵀ
164 M = SIZE2 (X); N = SIZE1 (X);
165 U = gsl_matrix_calloc (M, N);
166 gsl_matrix_transpose_memcpy (U, X);
167 } else {
168 // X = USVᵀ
169 U = gsl_matrix_calloc (M, N);
170 gsl_matrix_memcpy (U, X);
171 }
172 // X = USVᵀ by Jacobi, more precise than Golub-Reinsch.
173 // GSL routine yields V, not Vᵀ, U develops in place.
174 gsl_matrix *V = gsl_matrix_calloc (N, N);
175 (*Sv) = gsl_vector_calloc (N);
176 ASSERT_GSL (gsl_linalg_SV_decomp_jacobi (U, V, (*Sv)));
177 // Return singular values, if required.
178 gsl_matrix *eigens = gsl_matrix_calloc (M, N);
179 if (transpose) {
180 ASSERT_GSL (gsl_matrix_memcpy (eigens, U));
181 } else {
182 ASSERT_GSL (gsl_matrix_memcpy (eigens, V));
183 }
184 a68_matrix_free (U);
185 a68_matrix_free (V);
186 return eigens;
187 }
188
189 //! @brief PROC pseudo inv = ([, ] REAL, REAL) [, ] REAL
190
191 void genie_matrix_pinv_lim (NODE_T * p)
192 {
193 // Compute the Moore-Penrose pseudo inverse of a MxN matrix.
194 // G. Strang. "Linear Algebra and its applications", 2nd edition.
195 // Academic Press [1980].
196 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
197 A68_REAL lim;
198 POP_OBJECT (p, &lim, A68_REAL);
199 gsl_matrix *X = pop_matrix (p, A68_TRUE), *mpinv = NO_REAL_MATRIX;
200 compute_pseudo_inverse (p, &mpinv, X, VALUE (&lim));
201 push_matrix (p, mpinv);
202 a68_matrix_free (mpinv);
203 a68_matrix_free (X);
204 (void) gsl_set_error_handler (save_handler);
205 }
206
207 //! @brief OP PINV = ([, ] REAL) [, ] REAL
208
209 void genie_matrix_pinv (NODE_T * p)
210 {
211 PUSH_VALUE (p, SMALL_EIGEN, A68_REAL);
212 genie_matrix_pinv_lim (p);
213 }
214
215 //! @brief OP MEAN = ([, ] REAL) [, ] REAL
216
217 void genie_matrix_column_mean (NODE_T *p)
218 {
219 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
220 gsl_matrix *Z = pop_matrix (p, A68_TRUE);
221 unt M = SIZE1 (Z), N = SIZE2 (Z);
222 gsl_matrix *A = gsl_matrix_calloc (M, N);
223 for (int i = 0; i < N; i++) {
224 DOUBLE_T sum = 0;
225 for (int j = 0; j < M; j++) {
226 sum += gsl_matrix_get (Z, j, i);
227 }
228 DOUBLE_T mean = sum / M;
229 for (int j = 0; j < M; j++) {
230 gsl_matrix_set (A, j, i, mean);
231 }
232 }
233 push_matrix (p, A);
234 a68_matrix_free (A);
235 a68_matrix_free (Z);
236 (void) gsl_set_error_handler (save_handler);
237 }
238
239 //! @brief PROC left columns = ([, ] REAL, INT) [, ] REAL
240
241 void genie_left_columns (NODE_T *p)
242 {
243 A68_INT k;
244 POP_OBJECT (p, &k, A68_INT);
245 gsl_matrix *X = pop_matrix (p, A68_TRUE);
246 gsl_matrix *Y = left_columns (p, X, VALUE (&k));
247 push_matrix (p, Y);
248 a68_matrix_free (X);
249 a68_matrix_free (Y);
250 }
251
252 //! @brief PROC pca cv = ([, ] REAL, REF [] REAL) [, ] REAL
253
254 void genie_matrix_pca_cv (NODE_T * p)
255 {
256 // Principal component analysis of a MxN matrix from a covariance matrix.
257 // Forming the covariance matrix squares the condition number.
258 // Hence this routine might loose more significant digits than SVD.
259 // On the other hand, using PCA one looks for dominant eigenvectors.
260 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
261 A68_REF singulars;
262 POP_REF (p, &singulars);
263 gsl_matrix *X = pop_matrix (p, A68_TRUE);
264 MATH_RTE (p, X == NO_REAL_MATRIX, M_ROW_ROW_REAL, "empty data matrix");
265 // M samples, N features.
266 unt M = SIZE1 (X), N = SIZE2 (X);
267 // Keep X column-mean-centered.
268 gsl_matrix *C = column_mean (X);
269 // Covariance matrix, M samples: Cov = XᵀX.
270 M = MAX (M, N);
271 gsl_matrix *CV = NO_REAL_MATRIX;
272 a68_dgemm (FLIP, SELF, 1, C, C, 0, &CV);
273 // Compute and sort eigenvectors.
274 gsl_vector *Sv = gsl_vector_calloc (M);
275 gsl_matrix *eigens = gsl_matrix_calloc (M, M);
276 gsl_eigen_symmv_workspace *W = gsl_eigen_symmv_alloc (M);
277 ASSERT_GSL (gsl_eigen_symmv (CV, Sv, eigens, W));
278 ASSERT_GSL (gsl_eigen_symmv_sort (Sv, eigens, GSL_EIGEN_SORT_ABS_DESC));
279 // Yield results.
280 if (!IS_NIL (singulars)) {
281 *DEREF (A68_REF, &singulars) = vector_to_row (p, Sv);
282 }
283 push_matrix (p, eigens);
284 (void) gsl_set_error_handler (save_handler);
285 // Garbage collector
286 a68_matrix_free (eigens);
287 a68_matrix_free (X);
288 gsl_eigen_symmv_free (W);
289 a68_matrix_free (C);
290 a68_matrix_free (CV);
291 a68_vector_free (Sv);
292 }
293
294 //! @brief PROC pca svd = ([, ] REAL, REF [] REAL) [, ] REAL
295
296 void genie_matrix_pca_svd (NODE_T * p)
297 {
298 // Principal component analysis of a MxN matrix by Jacobi SVD.
299 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
300 A68_REF singulars;
301 POP_REF (p, &singulars);
302 // X data matrix, samples in rows, features in columns
303 gsl_matrix *X = pop_matrix (p, A68_TRUE);
304 MATH_RTE (p, X == NO_REAL_MATRIX, M_ROW_ROW_REAL, "empty data matrix");
305 // Keep X column-mean-centered.
306 gsl_matrix *C = column_mean (X);
307 gsl_vector *Sv = NO_REAL_VECTOR;
308 gsl_matrix *eigens = pca_svd (C, &Sv);
309 // Yield results.
310 if (!IS_NIL (singulars)) {
311 *DEREF (A68_REF, &singulars) = vector_to_row (p, Sv);
312 }
313 push_matrix (p, eigens);
314 (void) gsl_set_error_handler (save_handler);
315 // Clean up.
316 a68_matrix_free (eigens);
317 a68_matrix_free (X);
318 a68_matrix_free (C);
319 if (Sv != NO_REAL_VECTOR) {
320 a68_vector_free (Sv);
321 }
322 }
323
324 //! @brief PROC pcr = ([, ] REAL, [, ] REAL, REF [] REAL, INT, REAL) [, ] REAL
325
326 void genie_matrix_pcr (NODE_T * p)
327 {
328 // Principal Component Regression of a MxN matrix by Jacobi SVD.
329 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
330 // Pop arguments.
331 // 'lim' is minimum relative length to first eigenvector.
332 A68_REAL lim;
333 POP_OBJECT (p, &lim, A68_REAL);
334 // 'Nc' is maximum number of eigenvectors.
335 A68_INT Nc;
336 POP_OBJECT (p, &Nc, A68_INT);
337 A68_REF singulars;
338 POP_REF (p, &singulars);
339 // Y data vector, features in a single column
340 gsl_matrix *Y = pop_matrix (p, A68_TRUE);
341 MATH_RTE (p, Y == NO_REAL_MATRIX, M_ROW_ROW_REAL, "empty data matrix");
342 unt My = SIZE1 (Y), Ny = SIZE2 (Y);
343 // X data matrix, samples in rows, features in columns
344 gsl_matrix *X = pop_matrix (p, A68_TRUE);
345 MATH_RTE (p, X == NO_REAL_MATRIX, M_ROW_ROW_REAL, "empty data matrix");
346 unt Mx = SIZE1 (X), Nx = SIZE2 (X);
347 // Catch wrong calls.
348 MATH_RTE (p, My == 0 || Ny == 0, M_ROW_ROW_REAL, "invalid column vector F");
349 MATH_RTE (p, Mx == 0 || Nx == 0, M_ROW_ROW_REAL, "invalid data matrix E");
350 MATH_RTE (p, Mx != My, M_ROW_ROW_REAL, "rows in F must match columns in E");
351 MATH_RTE (p, Ny != 1, M_ROW_ROW_REAL, "F must be a column vector");
352 MATH_RTE (p, VALUE (&lim) < 0 || VALUE (&lim) > 1.0, M_REAL, "invalid relative length");
353 // Keep X column-mean-centered.
354 gsl_matrix *C = column_mean (X);
355 gsl_vector *Sv = NO_REAL_VECTOR;
356 gsl_matrix *eigens = pca_svd (C, &Sv);
357 // Dimension reduction.
358 int Nk = (VALUE (&Nc) == 0 ? SIZE (Sv) : MIN (SIZE (Sv), VALUE (&Nc)));
359 if (VALUE (&lim) <= 0) {
360 VALUE (&lim) = SMALL_EIGEN;
361 }
362 // Too short eigenvectors wreak havoc.
363 REAL_T Sv0 = gsl_vector_get (Sv, 0);
364 for (int k = 1; k < SIZE (Sv); k++) {
365 if (gsl_vector_get (Sv, k) / Sv0 < VALUE (&lim)) {
366 if (k + 1 < Nk) {
367 Nk = k + 1;
368 }
369 }
370 }
371 gsl_matrix *reduced = left_columns (p, eigens, Nk);
372 // Compute projected set = X * reduced.
373 gsl_matrix *proj = NO_REAL_MATRIX;
374 a68_dgemm (SELF, SELF, 1.0, X, reduced, 0.0, &proj);
375 // Compute β = reduced * P⁻¹ * Y.
376 gsl_matrix *mpinv = NO_REAL_MATRIX;
377 compute_pseudo_inverse (p, &mpinv, proj, SMALL_EIGEN);
378 gsl_matrix *z = NO_REAL_MATRIX, *beta = NO_REAL_MATRIX;
379 a68_dgemm (SELF, SELF, 1.0, reduced, mpinv, 0.0, &z);
380 a68_dgemm (SELF, SELF, 1.0, z, Y, 0.0, &beta);
381 // Yield results.
382 if (!IS_NIL (singulars)) {
383 *DEREF (A68_REF, &singulars) = vector_to_row (p, Sv);
384 }
385 push_matrix (p, beta);
386 (void) gsl_set_error_handler (save_handler);
387 // Clean up.
388 a68_matrix_free (eigens);
389 a68_matrix_free (reduced);
390 a68_matrix_free (X);
391 a68_matrix_free (Y);
392 a68_matrix_free (C);
393 a68_vector_free (Sv);
394 a68_matrix_free (z);
395 a68_matrix_free (mpinv);
396 a68_matrix_free (proj);
397 a68_matrix_free (beta);
398 }
399
400 //! @brief PROC ols = ([, ] REAL, [, ] REAL) [, ] REAL
401
402 void genie_matrix_ols (NODE_T *p)
403 {
404 // TLS relates to PLS as OLS relates to PCR.
405 // Note that X is an MxN matrix and Y, β are Nx1 column vectors.
406 // OLS can implemented (inefficiently) as PCR with maximum number of singular values:
407 // PUSH_VALUE (p, 0.0, A68_REAL);
408 // PUSH_VALUE (p, 0, A68_INT);
409 // genie_matrix_pcr (p);
410 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
411 // X and Y matrices.
412 gsl_matrix *Y = pop_matrix (p, A68_TRUE);
413 gsl_matrix *X = pop_matrix (p, A68_TRUE);
414 // Compute β = X⁻¹ * Y.
415 gsl_matrix *mpinv = NO_REAL_MATRIX, *beta = NO_REAL_MATRIX;
416 compute_pseudo_inverse (p, &mpinv, X, SMALL_EIGEN);
417 a68_dgemm (SELF, SELF, 1.0, X, Y, 0.0, &beta);
418 // Yield results.
419 push_matrix (p, beta);
420 (void) gsl_set_error_handler (save_handler);
421 a68_matrix_free (beta);
422 a68_matrix_free (mpinv);
423 a68_matrix_free (X);
424 a68_matrix_free (Y);
425 }
426
427 //! @brief PROC pls1 = ([, ] REAL, [, ] REAL, REF [] REAL, INT, REAL) [, ] REAL
428
429 void genie_matrix_pls1 (NODE_T *p)
430 {
431 // PLS decomposes X and Y data concurrently as
432 // X = T Pᵀ + E
433 // Y = U Qᵀ + F
434 // PLS1 is a widely used algorithm appropriate for the vector Y case.
435 // This procedure computes PLS1 with SVD solver for β,
436 // by a NIPALS algorithm with orthonormal scores and loadings.
437 // See Ulf Indahl, Journal of Chemometrics 28(3) 168-180 [2014].
438 // Note that E is an MxN matrix and F, β are Nx1 column vectors.
439 // This for consistency with PLS2.
440 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
441 // Pop arguments.
442 // 'lim' is minimum relative length to first eigenvector.
443 A68_REAL lim;
444 POP_OBJECT (p, &lim, A68_REAL);
445 // 'Nc' is maximum number of eigenvectors.
446 A68_INT Nc;
447 POP_OBJECT (p, &Nc, A68_INT);
448 A68_REF singulars;
449 POP_REF (p, &singulars);
450 // X and Y matrices.
451 gsl_matrix *F = pop_matrix (p, A68_TRUE);
452 gsl_matrix *E = pop_matrix (p, A68_TRUE);
453 // Catch wrong calls.
454 unt Me = SIZE1 (E), Ne = SIZE2 (E), Mf = SIZE1 (F), Nf = SIZE2 (F);
455 MATH_RTE (p, Mf == 0 || Nf == 0, M_ROW_ROW_REAL, "invalid column vector F");
456 MATH_RTE (p, Me == 0 || Ne == 0, M_ROW_ROW_REAL, "invalid data matrix E");
457 MATH_RTE (p, Me != Mf, M_ROW_ROW_REAL, "rows in F must match columns in E");
458 MATH_RTE (p, Nf != 1, M_ROW_ROW_REAL, "F must be a column vector");
459 MATH_RTE (p, VALUE (&lim) < 0 || VALUE (&lim) > 1.0, M_REAL, "invalid relative length");
460 // Sensible defaults.
461 int Nk = (VALUE (&Nc) == 0 ? MIN (Ne, Mf) : MIN (Mf, VALUE (&Nc)));
462 if (VALUE (&lim) <= 0) {
463 VALUE (&lim) = SMALL_EIGEN;
464 }
465 // Decompose E and F.
466 gsl_matrix *EIGEN = NO_REAL_MATRIX, *nE = NO_REAL_MATRIX, *nF = NO_REAL_MATRIX;
467 gsl_vector *Sv = gsl_vector_calloc (Nk);
468 BOOL_T siga = A68_TRUE;
469 gsl_matrix *eigens = NO_REAL_MATRIX, *lat = NO_REAL_MATRIX, *pl = NO_REAL_MATRIX, *ql = NO_REAL_MATRIX;
470 for (int k = 0; k < Nk && siga; k ++) {
471 // E weight from E, F covariance.
472 // eigens = Eᵀ * f / | Eᵀ * f |
473 a68_dgemm (FLIP, SELF, 1.0, E, F, 0.0, &eigens);
474 REAL_T norm = matrix_norm (eigens);
475 if (k > 0 && (norm / gsl_vector_get (Sv, 0)) < VALUE (&lim)) {
476 Nk = k;
477 siga = A68_FALSE;
478 } else {
479 ASSERT_GSL (gsl_matrix_scale (eigens, 1.0 / norm));
480 gsl_vector_set (Sv, k, norm);
481 // Compute latent variable.
482 // lat = E * eigens / | E * eigens |
483 a68_dgemm (SELF, SELF, 1.0, E, eigens, 0.0, &lat);
484 norm = matrix_norm (lat);
485 ASSERT_GSL (gsl_matrix_scale (lat, 1.0 / norm));
486 // Deflate E and F, remove latent variable from both.
487 // pl = Eᵀ * lat; E -= lat * plᵀ and ql = Fᵀ * lat; F -= lat * qlᵀ
488 a68_dgemm (FLIP, SELF, 1.0, E, lat, 0.0, &pl);
489 a68_dgemm (SELF, FLIP, -1.0, lat, pl, 1.0, &E);
490 a68_dgemm (FLIP, SELF, 1.0, F, lat, 0.0, &ql);
491 a68_dgemm (SELF, FLIP, -1.0, lat, ql, 1.0, &F);
492 // Build matrices.
493 EIGEN = mat_before_ab (p, EIGEN, eigens);
494 nE = mat_before_ab (p, nE, pl); // P
495 nF = mat_over_ab (p, nF, ql); // Qᵀ
496 }
497 }
498 // Projection of original data = Eᵀ * EIGEN
499 gsl_matrix *nP = NO_REAL_MATRIX;
500 a68_dgemm (FLIP, SELF, 1.0, nE, EIGEN, 0.0, &nP);
501 // Compute β = EIGEN * nP⁻¹ * nF
502 gsl_matrix *mpinv = NO_REAL_MATRIX, *z = NO_REAL_MATRIX, *beta = NO_REAL_MATRIX;
503 compute_pseudo_inverse (p, &mpinv, nP, SMALL_EIGEN);
504 a68_dgemm (SELF, SELF, 1.0, EIGEN, mpinv, 0.0, &z);
505 a68_dgemm (SELF, SELF, 1.0, z, nF, 0.0, &beta);
506 // Yield results.
507 if (!IS_NIL (singulars)) {
508 gsl_vector *Svl = gsl_vector_calloc (Nk);
509 for (int k = 0; k < Nk; k++) {
510 gsl_vector_set (Svl, k, gsl_vector_get (Sv, k));
511 }
512 *DEREF (A68_REF, &singulars) = vector_to_row (p, Svl);
513 a68_vector_free (Svl);
514 }
515 push_matrix (p, beta);
516 (void) gsl_set_error_handler (save_handler);
517 // Garbage collector.
518 a68_matrix_free (E);
519 a68_matrix_free (F);
520 a68_matrix_free (eigens);
521 a68_matrix_free (lat);
522 a68_matrix_free (pl);
523 a68_matrix_free (ql);
524 a68_matrix_free (beta);
525 a68_matrix_free (EIGEN);
526 a68_matrix_free (nE);
527 a68_matrix_free (nF);
528 a68_matrix_free (nP);
529 a68_vector_free (Sv);
530 }
531
532 //! @brief PROC pls2 = ([, ] REAL, [, ] REAL, REF [] REAL, INT, REAL) [, ] REAL
533
534 void genie_matrix_pls2 (NODE_T *p)
535 {
536 // PLS decomposes X and Y data concurrently as
537 // X = T Pᵀ + E
538 // Y = U Qᵀ + F
539 // Generic orthodox NIPALS, following Hervé Abdi, University of Texas.
540 // "Partial Least Squares Regression and Projection on Latent Structure Regression",
541 // Computational Statistics [2010].
542 #define PLS_MAX_ITER 100
543 gsl_error_handler_t *save_handler = gsl_set_error_handler (torrix_error_handler);
544 // Pop arguments.
545 // 'lim' is minimum relative length to first eigenvector.
546 A68_REAL lim;
547 POP_OBJECT (p, &lim, A68_REAL);
548 // 'Nc' is maximum number of eigenvectors.
549 A68_INT Nc;
550 POP_OBJECT (p, &Nc, A68_INT);
551 A68_REF singulars;
552 POP_REF (p, &singulars);
553 // X and Y matrices.
554 gsl_matrix *F = pop_matrix (p, A68_TRUE);
555 gsl_matrix *E = pop_matrix (p, A68_TRUE);
556 // Catch wrong calls.
557 unt Me = SIZE1 (E), Ne = SIZE2 (E), Mf = SIZE1 (F), Nf = SIZE2 (F);
558 MATH_RTE (p, Mf == 0 || Nf == 0, M_ROW_ROW_REAL, "invalid column vector F");
559 MATH_RTE (p, Me == 0 || Ne == 0, M_ROW_ROW_REAL, "invalid data matrix E");
560 MATH_RTE (p, Me != Mf, M_ROW_ROW_REAL, "rows in F must match columns in E");
561 CHECK_INT_SHORTEN(p, VALUE (&Nc));
562 // Sensible defaults.
563 int Nk = (VALUE (&Nc) == 0 ? MIN (Ne, Mf) : MIN (Mf, VALUE (&Nc)));
564 if (VALUE (&lim) <= 0) {
565 VALUE (&lim) = SMALL_EIGEN;
566 }
567 // Initial vector u.
568 gsl_matrix *u = gsl_matrix_calloc (Mf, 1);
569 if (Nf == 1) { // Column vector, PLS1
570 for (int k = 0; k < Mf; k ++) {
571 gsl_matrix_set (u, k, 0, gsl_matrix_get (F, k, 0));
572 }
573 } else {
574 for (int k = 0; k < Mf; k ++) {
575 gsl_matrix_set (u, k, 0, a68_gauss_rand ());
576 }
577 }
578 // Decompose E and F jointly.
579 gsl_vector *Sv = gsl_vector_calloc (Nk), *dD = gsl_vector_calloc (Nk);
580 gsl_matrix *u0 = gsl_matrix_calloc (Mf, 1), *diff = gsl_matrix_calloc (Mf, 1);
581 gsl_matrix *eigens = NO_REAL_MATRIX;
582 gsl_matrix *t = NO_REAL_MATRIX, *b = NO_REAL_MATRIX, *c = NO_REAL_MATRIX;
583 gsl_matrix *pl = NO_REAL_MATRIX, *ql = NO_REAL_MATRIX, *nP = NO_REAL_MATRIX, *nC = NO_REAL_MATRIX;
584 BOOL_T siga = A68_TRUE;
585 for (int k = 0; k < Nk && siga; k ++) {
586 REAL_T norm_e, norm;
587 for (int j = 0; j < PLS_MAX_ITER; j ++) {
588 gsl_matrix_memcpy (u0, u);
589 // Compute X weight. Note that Eᵀ * u is of form [Indahl]
590 // Eᵀ * F * Fᵀ * E * w / |Eᵀ * F * Fᵀ * w| = (1 / lambda) (Fᵀ * E)ᵀ * (Fᵀ * E) * w,
591 // that is, w is an eigenvector of covariance matrix Fᵀ * E and
592 // longest eigenvector of symmetric matrix Eᵀ * F * Fᵀ * E.
593 a68_dgemm (FLIP, SELF, 1.0, E, u, 0.0, &eigens);
594 norm_e = matrix_norm (eigens);
595 ASSERT_GSL (gsl_matrix_scale (eigens, 1.0 / norm_e));
596 // X factor score.
597 a68_dgemm (SELF, SELF, 1.0, E, eigens, 0.0, &t);
598 norm = matrix_norm (t);
599 ASSERT_GSL (gsl_matrix_scale (t, 1.0 / norm));
600 // Y weight.
601 a68_dgemm (FLIP, SELF, 1.0, F, t, 0.0, &c);
602 norm = matrix_norm (c);
603 ASSERT_GSL (gsl_matrix_scale (c, 1.0 / norm));
604 // Y score.
605 a68_dgemm (SELF, SELF, 1.0, F, c, 0.0, &u);
606 gsl_matrix_memcpy (diff, u);
607 gsl_matrix_sub (diff, u0);
608 norm = matrix_norm (diff);
609 if (norm < SMALL_EIGEN) {
610 j = PLS_MAX_ITER;
611 }
612 }
613 if (k > 0 && (norm_e / gsl_vector_get (Sv, 0)) < VALUE (&lim)) {
614 Nk = k;
615 siga = A68_FALSE;
616 } else {
617 gsl_vector_set (Sv, k, norm_e);
618 // X factor loading and deflation.
619 a68_dgemm (FLIP, SELF, 1.0, E, t, 0.0, &pl);
620 norm = matrix_norm (t);
621 ASSERT_GSL (gsl_matrix_scale (pl, 1.0 / (norm * norm)));
622 a68_dgemm (SELF, FLIP, -1.0, t, pl, 1.0, &E);
623 // Y factor loading and deflation.
624 a68_dgemm (FLIP, SELF, 1.0, F, u, 0.0, &ql);
625 norm = matrix_norm (u);
626 ASSERT_GSL (gsl_matrix_scale (ql, 1.0 / (norm * norm)));
627 a68_dgemm (FLIP, SELF, 1.0, t, u, 0.0, &b);
628 a68_dgemm (SELF, FLIP, -gsl_matrix_get (b, 0, 0), t, ql, 1.0, &F);
629 // Build vector and matrices.
630 gsl_vector_set (dD, k, gsl_matrix_get(b, 0, 0));
631 nP = mat_before_ab (p, nP, pl); // P
632 nC = mat_before_ab (p, nC, c); // C
633 }
634 }
635 // Compute β = (Pᵀ)⁻¹ * D * Cᵀ
636 // Python diag function.
637 gsl_matrix *D = gsl_matrix_calloc (Nk, Nk);
638 for (int k = 0; k < Nk; k++) {
639 gsl_matrix_set (D, k, k, gsl_vector_get (dD, k));
640 }
641 gsl_matrix *mpinv = NO_REAL_MATRIX;
642 compute_pseudo_inverse (p, &mpinv, nP, SMALL_EIGEN);
643 gsl_matrix *z = NO_REAL_MATRIX, *beta = NO_REAL_MATRIX;
644 a68_dgemm (FLIP, SELF, 1.0, mpinv, D, 0.0, &z);
645 a68_dgemm (SELF, FLIP, 1.0, z, nC, 0.0, &beta);
646 // Yield results.
647 if (!IS_NIL (singulars)) {
648 gsl_vector *Svl = gsl_vector_calloc (Nk);
649 for (int k = 0; k < Nk; k++) {
650 gsl_vector_set (Svl, k, gsl_vector_get (Sv, k));
651 }
652 *DEREF (A68_REF, &singulars) = vector_to_row (p, Svl);
653 a68_vector_free (Svl);
654 }
655 push_matrix (p, beta);
656 (void) gsl_set_error_handler (save_handler);
657 // Garbage collector.
658 a68_vector_free (dD);
659 a68_vector_free (Sv);
660 a68_matrix_free (b);
661 a68_matrix_free (beta);
662 a68_matrix_free (c);
663 a68_matrix_free (D);
664 a68_matrix_free (diff);
665 a68_matrix_free (eigens);
666 a68_matrix_free (nC);
667 a68_matrix_free (nP);
668 a68_matrix_free (mpinv);
669 a68_matrix_free (pl);
670 a68_matrix_free (ql);
671 a68_matrix_free (t);
672 a68_matrix_free (u);
673 a68_matrix_free (u0);
674 a68_matrix_free (z);
675 #undef PLS_MAX_ITER
676 }
677
678 //! @brief PROC tls = ([, ] REAL, [, ] REAL, REF [] REAL) [, ] REAL
679
680 void genie_matrix_tls (NODE_T *p)
681 {
682 // TLS relates to PLS as OLS relates to PCR.
683 // TLS is implemented as PLS1 with maximum number of singular values.
684 PUSH_VALUE (p, 0.0, A68_REAL);
685 PUSH_VALUE (p, 0, A68_INT);
686 genie_matrix_pls1 (p);
687 }
688
689 #endif
