MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
matrix_svd.h
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Daniel Thuerck, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 *
9 * The matrix formats for this implementation are exemplary visualized:
10 *
11 * A A U U
12 * A A = U U * S S * V V
13 * A A U U S S V V
14 *
15 * S S S V V V
16 * A A A = U U U * S S S * V V V
17 * A A A U U U S S S V V V
18 */
19
20#ifndef MATH_MATRIX_SVD_HEADER
21#define MATH_MATRIX_SVD_HEADER
22
23#include <vector>
24
25#include "math/defines.h"
26#include "math/matrix.h"
27#include "math/matrix_tools.h"
28#include "math/matrix_qr.h"
29
30#define MATH_SVD_DEFAULT_ZERO_THRESHOLD 1e-12
31
33
49template <typename T>
50void
51matrix_svd (T const* mat_a, int rows, int cols,
52 T* mat_u, T* vec_s, T* mat_v,
53 T const& epsilon = T(MATH_SVD_DEFAULT_ZERO_THRESHOLD));
54
60template <typename T, int M, int N>
61void
62matrix_svd (Matrix<T, M, N> const& mat_a, Matrix<T, M, N>* mat_u,
63 Matrix<T, N, N>* mat_s, Matrix<T, N, N>* mat_v,
64 T const& epsilon = T(MATH_SVD_DEFAULT_ZERO_THRESHOLD));
65
72template <typename T, int M, int N>
73void
74matrix_pseudo_inverse (Matrix<T, M, N> const& A,
75 Matrix<T, N, M>* result,
76 T const& epsilon = T(MATH_SVD_DEFAULT_ZERO_THRESHOLD));
77
79
80/* ------------------------ SVD Internals ------------------------- */
81
84
91template <typename T>
92bool
93matrix_is_submatrix_zero_enclosed (T const* mat, int m, int k, T const& epsilon)
94{
95 int const j = m - k - 1;
96 if (j < 0)
97 return true;
98 for (int i = m - k; i < m; ++i)
99 if (!MATH_EPSILON_EQ(mat[j * m + i], T(0), epsilon)
100 || !MATH_EPSILON_EQ(mat[i * m + j], T(0), epsilon))
101 return false;
102 return true;
103}
104
109template <typename T>
110bool
112 int rows, int cols, T const& epsilon)
113{
114 int const n = std::min(rows, cols) - 1;
115 for (int i = 0; i < n; ++i)
116 if (MATH_EPSILON_EQ(T(0), mat[i * cols + i + 1], epsilon))
117 return false;
118 return true;
119}
120
125template <typename T>
126void
127matrix_2x2_eigenvalues (T const* mat, T* smaller_ev, T* larger_ev)
128{
129 /* For matrix [a b; c d] solve (a+d) / 2 + sqrt((a+d)^2 / 4 - ad + bc). */
130 T const& a = mat[0];
131 T const& b = mat[1];
132 T const& c = mat[2];
133 T const& d = mat[3];
134
135 T x = MATH_POW2(a + d) / T(4) - a * d + b * c;
136 x = (x > T(0) ? std::sqrt(x) : T(0));
137 *smaller_ev = (a + d) / T(2) - x;
138 *larger_ev = (a + d) / T(2) + x;
139}
140
147template <typename T>
148void
149matrix_householder_vector (T const* input, int length,
150 T* vector, T* beta, T const& epsilon, T const& norm_factor)
151{
152 T sigma(0);
153 for (int i = 1; i < length; ++i)
154 sigma += MATH_POW2(input[i] / norm_factor);
155
156 vector[0] = T(1);
157 for (int i = 1; i < length; ++i)
158 vector[i] = input[i] / norm_factor;
159
160 if (MATH_EPSILON_EQ(sigma, T(0), epsilon))
161 {
162 *beta = T(0);
163 return;
164 }
165
166 T first = input[0] / norm_factor;
167 T mu = std::sqrt(MATH_POW2(first) + sigma);
168 if (first < epsilon)
169 vector[0] = first - mu;
170 else
171 vector[0] = -sigma / (first + mu);
172
173 first = vector[0];
174 *beta = T(2) * MATH_POW2(first) / (sigma + MATH_POW2(first));
175 for (int i = 0; i < length; ++i)
176 vector[i] /= first;
177}
178
184template <typename T>
185void
186matrix_householder_matrix (T const* vector, int length, T const beta,
187 T* matrix)
188{
189 std::fill(matrix, matrix + length * length, T(0));
190 for (int i = 0; i < length; ++i)
191 {
192 matrix[i * length + i] = T(1);
193 }
194
195 for (int i = 0; i < length; ++i)
196 {
197 for (int j = 0; j < length; ++j)
198 {
199 matrix[i * length + j] -= beta * vector[i] * vector[j];
200 }
201 }
202}
203
208template <typename T>
209void
210matrix_apply_householder_matrix (T* mat_a, int rows, int cols,
211 T const* house_mat, int house_length, int offset_rows, int offset_cols)
212{
213 /* Save block from old matrix that will be modified. */
214 int house_length_n = house_length - (rows - cols);
215 std::vector<T> rhs(house_length * house_length_n);
216 for (int i = 0; i < house_length; ++i)
217 {
218 for (int j = 0; j < house_length_n; ++j)
219 {
220 rhs[i * house_length_n + j] = mat_a[(offset_rows + i) *
221 cols + (offset_cols + j)];
222 }
223 }
224
225 /* Multiply block matrices. */
226 for (int i = 0; i < (rows - offset_rows); ++i)
227 {
228 for (int j = 0; j < (cols - offset_cols); ++j)
229 {
230 T current(0);
231 for (int k = 0; k < house_length; ++k)
232 {
233 current += (house_mat[i * house_length + k] *
234 rhs[k * house_length_n + j]);
235 }
236 mat_a[(offset_rows + i) * cols + (offset_cols + j)] = current;
237 }
238 }
239}
240
248template <typename T>
249void
250matrix_bidiagonalize (T const* mat_a, int rows, int cols, T* mat_u,
251 T* mat_b, T* mat_v, T const& epsilon)
252{
253 /* Initialize U and V with identity matrices. */
254 matrix_set_identity(mat_u, rows);
255 matrix_set_identity(mat_v, cols);
256
257 /* Copy mat_a into mat_b. */
258 std::copy(mat_a, mat_a + rows * cols, mat_b);
259
260 int const steps = (rows == cols) ? (cols - 1) : cols;
261 for (int k = 0; k < steps; ++k)
262 {
263 int const sub_length = rows - k;
264 std::vector<T> buffer(sub_length // input_vec
265 + sub_length // house_vec
266 + sub_length * sub_length // house_mat
267 + rows * rows // update_u
268 + rows * rows); // mat_u_tmp
269 T* input_vec = &buffer[0];
270 T* house_vec = input_vec + sub_length;
271 T* house_mat = house_vec + sub_length;
272 T house_beta;
273 T* update_u = house_mat + sub_length * sub_length;
274 T* mat_u_tmp = update_u + rows * rows;
275
276 for (int i = 0; i < sub_length; ++i)
277 input_vec[i] = mat_b[(k + i) * cols + k];
278
279 matrix_householder_vector(input_vec, sub_length, house_vec,
280 &house_beta, epsilon, T(1));
281 matrix_householder_matrix(house_vec, sub_length,
282 house_beta, house_mat);
283 matrix_apply_householder_matrix(mat_b, rows, cols,
284 house_mat, sub_length, k, k);
285
286 for (int i = k + 1; i < rows; ++i)
287 mat_b[i * cols + k] = T(0);
288
289 /* Construct U update matrix and update U. */
290 std::fill(update_u, update_u + rows * rows, T(0));
291 for (int i = 0; i < k; ++i)
292 update_u[i * rows + i] = T(1);
293 for (int i = 0; i < sub_length; ++i)
294 {
295 for (int j = 0; j < sub_length; ++j)
296 {
297 update_u[(k + i) * rows + (k + j)] =
298 house_mat[i * sub_length + j];
299 }
300 }
301
302 /* Copy matrix U for multiplication. */
303 std::copy(mat_u, mat_u + rows * rows, mat_u_tmp);
304 matrix_multiply(mat_u_tmp, rows, rows, update_u, rows, mat_u);
305
306 if (k <= cols - 3)
307 {
308 /* Normalization constant for numerical stability. */
309 T norm(0);
310 for (int i = k + 1; i < cols; ++i)
311 norm += mat_b[k * cols + i];
312 if (MATH_EPSILON_EQ(norm, T(0), epsilon))
313 norm = T(1);
314
315 int const inner_sub_length = cols - (k + 1);
316 int const slice_rows = rows - k;
317 int const slice_cols = cols - k - 1;
318 std::vector<T> buffer2(inner_sub_length // inner_input_vec
319 + inner_sub_length // inner_house_vec
320 + inner_sub_length * inner_sub_length // inner_house_mat
321 + slice_rows * slice_cols // mat_b_res
322 + slice_rows * slice_cols // mat_b_tmp
323 + cols * cols // update_v
324 + cols * cols); // mat_v_tmp
325
326 T* inner_input_vec = &buffer2[0];
327 T* inner_house_vec = inner_input_vec + inner_sub_length;
328 T* inner_house_mat = inner_house_vec + inner_sub_length;
329 T inner_house_beta;
330 T* mat_b_res = inner_house_mat + inner_sub_length * inner_sub_length;
331 T* mat_b_tmp = mat_b_res + slice_rows * slice_cols;
332 T* update_v = mat_b_tmp + slice_rows * slice_cols;
333 T* mat_v_tmp = update_v + cols * cols;
334
335 for (int i = 0; i < cols - k - 1; ++i)
336 inner_input_vec[i] = mat_b[k * cols + (k + 1 + i)];
337
338 matrix_householder_vector(inner_input_vec, inner_sub_length,
339 inner_house_vec, &inner_house_beta, epsilon, norm);
340 matrix_householder_matrix(inner_house_vec, inner_sub_length,
341 inner_house_beta, inner_house_mat);
342
343 /* Cut out mat_b(k:m, (k+1):n). */
344 for (int i = 0; i < slice_rows; ++i)
345 {
346 for (int j = 0; j < slice_cols; ++j)
347 {
348 mat_b_tmp[i * slice_cols + j] =
349 mat_b[(k + i) * cols + (k + 1 + j)];
350 }
351 }
352 matrix_multiply(mat_b_tmp, slice_rows, slice_cols,
353 inner_house_mat, inner_sub_length, mat_b_res);
354
355 /* Write frame back into mat_b. */
356 for (int i = 0; i < slice_rows; ++i)
357 {
358 for (int j = 0; j < slice_cols; ++j)
359 {
360 mat_b[(k + i) * cols + (k + 1 + j)] =
361 mat_b_res[i * slice_cols + j];
362 }
363 }
364
365 for (int i = k + 2; i < cols; ++i)
366 mat_b[k * cols + i] = T(0);
367
368 std::fill(update_v, update_v + cols * cols, T(0));
369 for (int i = 0; i < k + 1; ++i)
370 update_v[i * cols + i] = T(1);
371 for (int i = 0; i < inner_sub_length; ++i)
372 {
373 for (int j = 0; j < inner_sub_length; ++j)
374 {
375 update_v[(k + i + 1) * cols + (k + j + 1)] =
376 inner_house_mat[i * inner_sub_length + j];
377 }
378 }
379
380 /* Copy matrix v for multiplication. */
381 std::copy(mat_v, mat_v + cols * cols, mat_v_tmp);
382 matrix_multiply(mat_v_tmp, cols, cols, update_v, cols, mat_v);
383 }
384 }
385}
386
390template <typename T>
391void
392matrix_gk_svd_step (int rows, int cols, T* mat_b, T* mat_q, T* mat_p,
393 int p, int q, T const& epsilon)
394{
395 int const slice_length = cols - q - p;
396 int const mat_sizes = slice_length * slice_length;
397 std::vector<T> buffer(3 * mat_sizes);
398 T* mat_b22 = &buffer[0];
399 T* mat_b22_t = mat_b22 + mat_sizes;
400 T* mat_tmp = mat_b22_t + mat_sizes;
401
402 for (int i = 0; i < slice_length; ++i)
403 for (int j = 0; j < slice_length; ++j)
404 mat_b22[i * slice_length + j] = mat_b[(p + i) * cols + (p + j)];
405 for (int i = 0; i < slice_length; ++i)
406 for (int j = 0; j < slice_length; ++j)
407 mat_b22_t[i * slice_length + j] = mat_b22[j * slice_length + i];
408
409 /* Slice outer product gives covariance matrix. */
410 matrix_multiply(mat_b22, slice_length, slice_length, mat_b22_t,
411 slice_length, mat_tmp);
412
413 T mat_c[2 * 2];
414 mat_c[0] = mat_tmp[(slice_length - 2) * slice_length + (slice_length - 2)];
415 mat_c[1] = mat_tmp[(slice_length - 2) * slice_length + (slice_length - 1)];
416 mat_c[2] = mat_tmp[(slice_length - 1) * slice_length + (slice_length - 2)];
417 mat_c[3] = mat_tmp[(slice_length - 1) * slice_length + (slice_length - 1)];
418
419 /* Use eigenvalue that is closer to the lower right entry of the slice. */
420 T eig_1, eig_2;
421 matrix_2x2_eigenvalues(mat_c, &eig_1, &eig_2);
422
423 T diff1 = std::abs(mat_c[3] - eig_1);
424 T diff2 = std::abs(mat_c[3] - eig_2);
425 T mu = (diff1 < diff2) ? eig_1 : eig_2;
426
427 /* Zero another entry bz applyting givens rotations. */
428 int k = p;
429 T alpha = mat_b[k * cols + k] * mat_b[k * cols + k] - mu;
430 T beta = mat_b[k * cols + k] * mat_b[k * cols + (k + 1)];
431
432 for (int k = p; k < cols - q - 1; ++k)
433 {
434 T givens_c, givens_s;
435 matrix_givens_rotation(alpha, beta, &givens_c, &givens_s, epsilon);
436 matrix_apply_givens_column(mat_b, cols, cols, k, k + 1, givens_c,
437 givens_s);
438 matrix_apply_givens_column(mat_p, cols, cols, k, k + 1, givens_c,
439 givens_s);
440
441 alpha = mat_b[k * cols + k];
442 beta = mat_b[(k + 1) * cols + k];
443 matrix_givens_rotation(alpha, beta, &givens_c, &givens_s, epsilon);
444 internal::matrix_apply_givens_row(mat_b, cols, cols, k, k + 1,
445 givens_c, givens_s);
446 internal::matrix_apply_givens_column(mat_q, rows, cols, k, k + 1,
447 givens_c, givens_s);
448
449 if (k < (cols - q - 2))
450 {
451 alpha = mat_b[k * cols + (k + 1)];
452 beta = mat_b[k * cols + (k + 2)];
453 }
454 }
455}
456
457template <typename T>
458void
459matrix_svd_clear_super_entry(int rows, int cols, T* mat_b, T* mat_q,
460 int row_index, T const& epsilon)
461{
462 for (int i = row_index + 1; i < cols; ++i)
463 {
464 if (MATH_EPSILON_EQ(mat_b[row_index * cols + i], T(0), epsilon))
465 {
466 mat_b[row_index * cols + i] = T(0);
467 break;
468 }
469
470 T norm = MATH_POW2(mat_b[row_index * cols + i])
471 + MATH_POW2(mat_b[i * cols + i]);
472 norm = std::sqrt(norm) * MATH_SIGN(mat_b[i * cols + i]);
473
474 T givens_c = mat_b[i * cols + i] / norm;
475 T givens_s = mat_b[row_index * cols + i] / norm;
476 matrix_apply_givens_row(mat_b, cols, cols, row_index,
477 i, givens_c, givens_s);
478 matrix_apply_givens_column(mat_q, rows, cols, row_index,
479 i, givens_c, givens_s);
480 }
481}
482
486template <typename T>
487void
488matrix_gk_svd (T const* mat_a, int rows, int cols,
489 T* mat_u, T* vec_s, T* mat_v, T const& epsilon)
490{
491 /* Allocate memory for temp matrices. */
492 int const mat_q_full_size = rows * rows;
493 int const mat_b_full_size = rows * cols;
494 int const mat_p_size = cols * cols;
495 int const mat_q_size = rows * cols;
496 int const mat_b_size = cols * cols;
497
498 std::vector<T> buffer(mat_q_full_size + mat_b_full_size
499 + mat_p_size + mat_q_size + mat_b_size);
500 T* mat_q_full = &buffer[0];
501 T* mat_b_full = mat_q_full + mat_q_full_size;
502 T* mat_p = mat_b_full + mat_b_full_size;
503 T* mat_q = mat_p + mat_p_size;
504 T* mat_b = mat_q + mat_q_size;
505
506 matrix_bidiagonalize(mat_a, rows, cols,
507 mat_q_full, mat_b_full, mat_p, epsilon);
508
509 /* Extract smaller matrices. */
510 for (int i = 0; i < rows; ++i)
511 {
512 for (int j = 0; j < cols; ++j)
513 {
514 mat_q[i * cols + j] = mat_q_full[i * rows + j];
515 }
516 }
517 std::copy(mat_b_full, mat_b_full + cols * cols, mat_b);
518
519 /* Avoid infinite loops and exit after maximum number of iterations. */
520 int const max_iterations = rows * cols;
521 int iteration = 0;
522 while (iteration < max_iterations)
523 {
524 iteration += 1;
525
526 /* Enforce exact zeros for numerical stability. */
527 for (int i = 0; i < cols * cols; ++i)
528 {
529 T const& entry = mat_b[i];
530 if (MATH_EPSILON_EQ(entry, T(0), epsilon))
531 mat_b[i] = T(0);
532 }
533
534 /* GK 2a. */
535 for (int i = 0; i < (cols - 1); ++i)
536 {
537 if (std::abs(mat_b[i * cols + (i + 1)]) <= epsilon *
538 std::abs(mat_b[i * cols + i] + mat_b[(i + 1) * cols + (i + 1)]))
539 {
540 mat_b[i * cols + (i + 1)] = T(0);
541 }
542 }
543
544 /* GK 2b. */
545 /* Select q such that b33 is diagonal and blocked by zeros. */
546 int q = 0;
547 for (int k = 0; k < cols; ++k)
548 {
549 int const slice_len = k + 1;
550 std::vector<T> mat_b33(slice_len * slice_len);
551 for (int i = 0; i < slice_len; ++i)
552 {
553 for (int j = 0; j < slice_len; ++j)
554 {
555 mat_b33[i * slice_len + j]
556 = mat_b[(cols - k - 1 + i) * cols + (cols - k - 1 + j)];
557 }
558 }
559
560 if (matrix_is_diagonal(&mat_b33[0], slice_len, slice_len, epsilon))
561 {
562 if (k < cols - 1)
563 {
565 mat_b, cols, k + 1, epsilon))
566 {
567 q = k + 1;
568 }
569 }
570 else
571 {
572 q = k + 1;
573 }
574 }
575 }
576
577 /* Select z := n-p-q such that B22 has no zero superdiagonal entry. */
578 int z = 0;
579 T* mat_b22_tmp = new T[(cols - q) * (cols - q)];
580 for (int k = 0; k < (cols - q); ++k)
581 {
582 int const slice_len = k + 1;
583 for (int i = 0; i < slice_len; ++i)
584 {
585 for (int j = 0; j < slice_len; ++j)
586 {
587 mat_b22_tmp[i * slice_len + j]
588 = mat_b[(cols - q - k - 1 + i)
589 * cols + (cols - q - k - 1 + j)];
590 }
591 }
593 mat_b22_tmp, slice_len, slice_len, epsilon))
594 {
595 z = k + 1;
596 }
597 }
598 delete[] mat_b22_tmp;
599
600 int const p = cols - q - z;
601
602 /* GK 2c. */
603 if (q == cols)
604 break;
605
606 bool diagonal_non_zero = true;
607 int nz = 0;
608 for (nz = p; nz < (cols - q - 1); ++nz)
609 {
610 if (MATH_EPSILON_EQ(mat_b[nz * cols + nz], T(0), epsilon))
611 {
612 diagonal_non_zero = false;
613 mat_b[nz * cols + nz] = T(0);
614 break;
615 }
616 }
617
618 if (diagonal_non_zero)
619 matrix_gk_svd_step(rows, cols, mat_b, mat_q, mat_p, p, q, epsilon);
620 else
621 matrix_svd_clear_super_entry(rows, cols, mat_b, mat_q, nz, epsilon);
622 }
623
624 /* Create resulting matrices and vector from temporary entities. */
625 std::copy(mat_q, mat_q + rows * cols, mat_u);
626 std::copy(mat_p, mat_p + cols * cols, mat_v);
627 for (int i = 0; i < cols; ++i)
628 vec_s[i] = mat_b[i * cols + i];
629
630 /* Correct signs. */
631 for (int i = 0; i < cols; ++i)
632 {
633 if (vec_s[i] < epsilon)
634 {
635 vec_s[i] = -vec_s[i];
636 for (int j = 0; j < rows; ++j)
637 {
638 int index = j * cols + i;
639 mat_u[index] = -mat_u[index];
640 }
641 }
642 }
643}
644
649template <typename T>
650void
651matrix_r_svd (T const* mat_a, int rows, int cols,
652 T* mat_u, T* vec_s, T* mat_v, T const& epsilon)
653{
654 /* Allocate memory for temp matrices. */
655 int const mat_q_size = rows * rows;
656 int const mat_r_size = rows * cols;
657 int const mat_u_tmp_size = rows * cols;
658 std::vector<T> buffer(mat_q_size + mat_r_size + mat_u_tmp_size);
659 T* mat_q = &buffer[0];
660 T* mat_r = mat_q + mat_q_size;
661 T* mat_u_tmp = mat_r + mat_r_size;
662
663 matrix_qr(mat_a, rows, cols, mat_q, mat_r, epsilon);
664 matrix_gk_svd(mat_r, cols, cols, mat_u_tmp, vec_s, mat_v, epsilon);
665 std::fill(mat_u_tmp + cols * cols, mat_u_tmp + rows * cols, T(0));
666
667 /* Adapt U for big matrices. */
668 matrix_multiply(mat_q, rows, rows, mat_u_tmp, cols, mat_u);
669}
670
675template <typename T>
676int
677find_largest_ev_index (T const* values, int length)
678{
679 T largest = T(0);
680 int index = -1;
681 for (int i = 0; i < length; ++i)
682 if (values[i] > largest)
683 {
684 largest = values[i];
685 index = i;
686 }
687 return index;
688}
689
692
693/* ------------------------ Implementation ------------------------ */
694
696
697template <typename T>
698void
699matrix_svd (T const* mat_a, int rows, int cols,
700 T* mat_u, T* vec_s, T* mat_v, T const& epsilon)
701{
702 /* Allow for null result matrices. */
703 std::vector<T> mat_u_tmp;
704 std::vector<T> vec_s_tmp;
705 if (vec_s == nullptr)
706 {
707 vec_s_tmp.resize(cols);
708 vec_s = &vec_s_tmp[0];
709 }
710 std::vector<T> mat_v_tmp;
711 if (mat_v == nullptr)
712 {
713 mat_v_tmp.resize(cols * cols);
714 mat_v = &mat_v_tmp[0];
715 }
716
717 /*
718 * Handle two cases: The regular case, where M >= N (rows >= cols), and
719 * the irregular case, where M < N (rows < cols). In the latter one,
720 * zero rows are appended to A until A is a square matrix.
721 */
722 if (rows >= cols)
723 {
724 /* Allow for null result U matrix. */
725 if (mat_u == nullptr)
726 {
727 mat_u_tmp.resize(rows * cols);
728 mat_u = &mat_u_tmp[0];
729 }
730
731 /* Perform economy SVD if rows > 5/3 cols to save some operations. */
732 if (rows >= 5 * cols / 3)
733 {
734 internal::matrix_r_svd(mat_a, rows, cols,
735 mat_u, vec_s, mat_v, epsilon);
736 }
737 else
738 {
739 internal::matrix_gk_svd(mat_a, rows, cols,
740 mat_u, vec_s, mat_v, epsilon);
741 }
742 }
743 else
744 {
745 /* Append zero rows to A until A is square. */
746 std::vector<T> mat_a_tmp(cols * cols, T(0));
747 std::copy(mat_a, mat_a + cols * rows, &mat_a_tmp[0]);
748
749 /* Temporarily resize U, allow for null result matrices. */
750 mat_u_tmp.resize(cols * cols);
751 internal::matrix_gk_svd(&mat_a_tmp[0], cols, cols,
752 &mat_u_tmp[0], vec_s, mat_v, epsilon);
753
754 /* Copy the result back to U leaving out the last rows. */
755 if (mat_u != nullptr)
756 std::copy(&mat_u_tmp[0], &mat_u_tmp[0] + rows * cols, mat_u);
757 else
758 mat_u = &mat_u_tmp[0];
759 }
760
761 /* Sort the eigenvalues in S and adapt the columns of U and V. */
762 for (int i = 0; i < cols; ++i)
763 {
764 int pos = internal::find_largest_ev_index(vec_s + i, cols - i);
765 if (pos < 0)
766 break;
767 if (pos == 0)
768 continue;
769 std::swap(vec_s[i], vec_s[i + pos]);
770 matrix_swap_columns(mat_u, rows, cols, i, i + pos);
771 matrix_swap_columns(mat_v, cols, cols, i, i + pos);
772 }
773}
774
775#if 0
776template <typename T>
777void
778matrix_svd (T const* mat_a, int rows, int cols,
779 T* mat_u, T* vec_s, T* mat_v, T const& epsilon)
780{
781 /* Allow for null result matrices. */
782 std::vector<T> mat_u_tmp;
783 std::vector<T> vec_s_tmp;
784 std::vector<T> mat_v_tmp;
785 if (mat_u == nullptr)
786 {
787 mat_u_tmp.resize(rows * cols);
788 mat_u = &mat_u_tmp[0];
789 }
790 if (vec_s == nullptr)
791 {
792 vec_s_tmp.resize(cols);
793 vec_s = &vec_s_tmp[0];
794 }
795 if (mat_v == nullptr)
796 {
797 mat_v_tmp.resize(cols * cols);
798 mat_v = &mat_v_tmp[0];
799 }
800
801 /*
802 * The SVD can only handle systems where rows > cols. Thus, in case
803 * of cols > rows, the input matrix is transposed and the resulting
804 * solution converted to the desired solution:
805 *
806 * A* = (U S V*)* = V S* U* = V S U*
807 *
808 * In order to output same-sized matrices for every problem, many zero
809 * rows and columns can appear in systems with cols > rows.
810 */
811 std::vector<T> mat_a_tmp;
812 bool was_transposed = false;
813 if (cols > rows)
814 {
815 mat_a_tmp.resize(rows * cols);
816 std::copy(mat_a, mat_a + rows * cols, mat_a_tmp.begin());
817 matrix_transpose(&mat_a_tmp[0], rows, cols);
818 std::swap(rows, cols);
819 std::swap(mat_u, mat_v);
820 mat_a = &mat_a_tmp[0];
821 was_transposed = true;
822 }
823
824 /* Perform economy SVD if rows > 5/3 cols to save some operations. */
825 if (rows >= 5 * cols / 3)
826 {
827 internal::matrix_r_svd(mat_a, rows, cols,
828 mat_u, vec_s, mat_v, epsilon);
829 }
830 else
831 {
832 internal::matrix_gk_svd(mat_a, rows, cols,
833 mat_u, vec_s, mat_v, epsilon);
834 }
835
836 if (was_transposed)
837 {
838 std::swap(rows, cols);
839 std::swap(mat_u, mat_v);
840
841 /* Fix S by appending zeros. */
842 for (int i = rows; i < cols; ++i)
843 vec_s[i] = T(0);
844
845 /* Fix U by reshaping the matrix. */
846 for (int i = rows * rows - 1, y = rows - 1; y >= 0; --y)
847 {
848 for (int x = rows - 1; x >= 0; --x, --i)
849 mat_u[y * cols + x] = mat_u[i];
850 for (int x = rows; x < cols; ++x)
851 mat_u[y * cols + x] = T(0);
852 }
853
854 /* Fix V by reshaping the matrix. */
855 for (int i = cols * rows - 1, y = cols - 1; y >= 0; --y)
856 {
857 for (int x = rows - 1; x >= 0; --x, --i)
858 mat_v[y * cols + x] = mat_v[i];
859 for (int x = rows; x < cols; ++x)
860 mat_v[y * cols + x] = T(0);
861 }
862 }
863
864 /* Sort the eigenvalues in S and adapt the columns of U and V. */
865 for (int i = 0; i < cols; ++i)
866 {
867 int pos = internal::find_largest_ev_index(vec_s + i, cols - i);
868 if (pos < 0)
869 break;
870 if (pos == 0)
871 continue;
872 std::swap(vec_s[i], vec_s[i + pos]);
873 matrix_swap_columns(mat_u, rows, cols, i, i + pos);
874 matrix_swap_columns(mat_v, cols, cols, i, i + pos);
875 }
876}
877#endif
878
879template <typename T, int M, int N>
880void
882 Matrix<T, N, N>* mat_s, Matrix<T, N, N>* mat_v, T const& epsilon)
883{
884 T* mat_u_ptr = mat_u ? mat_u->begin() : nullptr;
885 T* mat_s_ptr = mat_s ? mat_s->begin() : nullptr;
886 T* mat_v_ptr = mat_v ? mat_v->begin() : nullptr;
887
888 matrix_svd<T>(mat_a.begin(), M, N,
889 mat_u_ptr, mat_s_ptr, mat_v_ptr, epsilon);
890
891 /* Swap elements of S into place. */
892 if (mat_s_ptr)
893 {
894 std::fill(mat_s_ptr + N, mat_s_ptr + N * N, T(0));
895 for (int x = 1, i = N + 1; x < N; ++x, i += N + 1)
896 std::swap(mat_s_ptr[x], mat_s_ptr[i]);
897 }
898}
899
900template <typename T, int M, int N>
901void
903 Matrix<T, N, M>* result, T const& epsilon)
904{
908 matrix_svd(A, &U, &S, &V, epsilon);
909
910 /* Invert diagonal of S. */
911 for (int i = 0; i < N; ++i)
912 if (MATH_EPSILON_EQ(S(i, i), T(0), epsilon))
913 S(i, i) = T(0);
914 else
915 S(i, i) = T(1) / S(i, i);
916
917 *result = V * S * U.transposed();
918}
919
921
922#endif /* MATH_MATRIX_SVD_HEADER */
Matrix class for arbitrary dimensions and types.
Definition matrix.h:54
Matrix< T, M, N > transposed(void) const
Returns a transposed copy of self by treating rows as columns.
Definition matrix.h:448
T * begin(void)
Definition matrix.h:506
#define MATH_INTERNAL_NAMESPACE_BEGIN
Definition defines.h:24
#define MATH_EPSILON_EQ(x, v, eps)
Definition defines.h:96
#define MATH_NAMESPACE_BEGIN
Definition defines.h:15
#define MATH_NAMESPACE_END
Definition defines.h:16
#define MATH_INTERNAL_NAMESPACE_END
Definition defines.h:25
#define MATH_SIGN(x)
Definition defines.h:93
#define MATH_POW2(x)
Definition defines.h:68
#define MATH_SVD_DEFAULT_ZERO_THRESHOLD
Definition matrix_svd.h:30
std::size_t first
Definition mesh_info.cc:46
bool matrix_is_submatrix_zero_enclosed(T const *mat, int m, int k, T const &epsilon)
Checks whether the lower-right square sub-matrix of size KxK is enclosed by zeros (up to some epsilon...
Definition matrix_svd.h:93
void matrix_gk_svd_step(int rows, int cols, T *mat_b, T *mat_q, T *mat_p, int p, int q, T const &epsilon)
Single step in the [GK-SVD] method.
Definition matrix_svd.h:392
void matrix_r_svd(T const *mat_a, int rows, int cols, T *mat_u, T *vec_s, T *mat_v, T const &epsilon)
Implementation of the [R-SVD] method, uses [GK-SVD] as solver for the reduced problem.
Definition matrix_svd.h:651
void matrix_gk_svd(T const *mat_a, int rows, int cols, T *mat_u, T *vec_s, T *mat_v, T const &epsilon)
Implementation of the [GK-SVD] method.
Definition matrix_svd.h:488
void matrix_2x2_eigenvalues(T const *mat, T *smaller_ev, T *larger_ev)
Returns the larger eigenvalue of the given 2x2 matrix.
Definition matrix_svd.h:127
void matrix_apply_householder_matrix(T *mat_a, int rows, int cols, T const *house_mat, int house_length, int offset_rows, int offset_cols)
Applies a given householder matrix to a frame in a given matrix with offset (offset_rows,...
Definition matrix_svd.h:210
int find_largest_ev_index(T const *values, int length)
Returns the index of the largest eigenvalue.
Definition matrix_svd.h:677
bool matrix_is_superdiagonal_nonzero(T const *mat, int rows, int cols, T const &epsilon)
Checks whether the super-diagonal (above the diagonal) of a MxN matrix does not contain zeros up to s...
Definition matrix_svd.h:111
void matrix_apply_givens_row(T *mat, int, int cols, int givens_i, int givens_k, T const &givens_c, T const &givens_s)
Applies a transposed Givens rotation for rows (givens_i, givens_k) by only rotating the required set ...
Definition matrix_qr.h:104
void matrix_givens_rotation(T const &alpha, T const &beta, T *givens_c, T *givens_s, T const &epsilon)
Calculates the Givens rotation coefficients c and s by solving [alpha beta] [c s;-c s] = [sqrt(alpha^...
Definition matrix_qr.h:56
void matrix_bidiagonalize(T const *mat_a, int rows, int cols, T *mat_u, T *mat_b, T *mat_v, T const &epsilon)
Bidiagonalizes a given MxN matrix, resulting in a MxN matrix U, a bidiagonal MxN matrix B and a NxN m...
Definition matrix_svd.h:250
void matrix_svd_clear_super_entry(int rows, int cols, T *mat_b, T *mat_q, int row_index, T const &epsilon)
Definition matrix_svd.h:459
void matrix_apply_givens_column(T *mat, int rows, int cols, int givens_i, int givens_k, T const &givens_c, T const &givens_s)
Applies a Givens rotation for columns (givens_i, givens_k) by only rotating the required set of colum...
Definition matrix_qr.h:86
void matrix_householder_vector(T const *input, int length, T *vector, T *beta, T const &epsilon, T const &norm_factor)
Creates a householder transformation vector and the coefficient for householder matrix creation.
Definition matrix_svd.h:149
void matrix_householder_matrix(T const *vector, int length, T const beta, T *matrix)
Given a Householder vector and beta coefficient, this function creates a transformation matrix to app...
Definition matrix_svd.h:186
void matrix_multiply(T const *mat_a, int rows_a, int cols_a, T const *mat_b, int cols_b, T *mat_res)
Matrix multiplication of dynamically sized dense matrices.
void matrix_pseudo_inverse(Matrix< T, M, N > const &A, Matrix< T, N, M > *result, T const &epsilon=T(1e-12))
Computes the Moore–Penrose pseudoinverse of matrix A using the SVD.
Definition matrix_svd.h:902
void matrix_qr(T const *mat_a, int rows, int cols, T *mat_q, T *mat_r, T const &epsilon=T(1e-12))
Calculates a QR decomposition for a given matrix A.
Definition matrix_qr.h:125
void matrix_swap_columns(T *const mat, int rows, int cols, int c1, int c2)
Swaps the columns c1 and c2 of matrix mat with dimension rows, cols.
void matrix_svd(T const *mat_a, int rows, int cols, T *mat_u, T *vec_s, T *mat_v, T const &epsilon=T(1e-12))
SVD for dynamic-size matrices A of size MxN (M rows, N columns).
Definition matrix_svd.h:699
Matrix< T, N, N > & matrix_set_identity(Matrix< T, N, N > *mat)
Sets the given square matrix to the identity matrix.
void matrix_transpose(T const *mat, int rows, int cols)
In-place transpose of a dynamically sized dense matrix.
bool matrix_is_diagonal(T *const mat, int rows, int cols, T const &epsilon=T(0))
Checks whether the input matrix is a diagonal matrix.
void swap(mve::Image< T > &a, mve::Image< T > &b)
Specialization of std::swap for efficient image swapping.
Definition image.h:478