20#ifndef MATH_MATRIX_SVD_HEADER
21#define MATH_MATRIX_SVD_HEADER
30#define MATH_SVD_DEFAULT_ZERO_THRESHOLD 1e-12
51matrix_svd (T
const* mat_a,
int rows,
int cols,
52 T* mat_u, T* vec_s, T* mat_v,
60template <
typename T,
int M,
int N>
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,
72template <
typename T,
int M,
int N>
74matrix_pseudo_inverse (Matrix<T, M, N>
const& A,
75 Matrix<T, N, M>* result,
95 int const j = m - k - 1;
98 for (
int i = m - k; i < m; ++i)
112 int rows,
int cols, T
const& epsilon)
114 int const n = std::min(rows, cols) - 1;
115 for (
int i = 0; i < n; ++i)
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;
150 T* vector, T* beta, T
const& epsilon, T
const& norm_factor)
153 for (
int i = 1; i < length; ++i)
154 sigma +=
MATH_POW2(input[i] / norm_factor);
157 for (
int i = 1; i < length; ++i)
158 vector[i] = input[i] / norm_factor;
166 T
first = input[0] / norm_factor;
169 vector[0] =
first - mu;
171 vector[0] = -sigma / (
first + mu);
175 for (
int i = 0; i < length; ++i)
189 std::fill(matrix, matrix + length * length, T(0));
190 for (
int i = 0; i < length; ++i)
192 matrix[i * length + i] = T(1);
195 for (
int i = 0; i < length; ++i)
197 for (
int j = 0; j < length; ++j)
199 matrix[i * length + j] -= beta * vector[i] * vector[j];
211 T
const* house_mat,
int house_length,
int offset_rows,
int offset_cols)
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)
218 for (
int j = 0; j < house_length_n; ++j)
220 rhs[i * house_length_n + j] = mat_a[(offset_rows + i) *
221 cols + (offset_cols + j)];
226 for (
int i = 0; i < (rows - offset_rows); ++i)
228 for (
int j = 0; j < (cols - offset_cols); ++j)
231 for (
int k = 0; k < house_length; ++k)
233 current += (house_mat[i * house_length + k] *
234 rhs[k * house_length_n + j]);
236 mat_a[(offset_rows + i) * cols + (offset_cols + j)] = current;
251 T* mat_b, T* mat_v, T
const& epsilon)
258 std::copy(mat_a, mat_a + rows * cols, mat_b);
260 int const steps = (rows == cols) ? (cols - 1) : cols;
261 for (
int k = 0; k < steps; ++k)
263 int const sub_length = rows - k;
264 std::vector<T> buffer(sub_length
266 + sub_length * sub_length
269 T* input_vec = &buffer[0];
270 T* house_vec = input_vec + sub_length;
271 T* house_mat = house_vec + sub_length;
273 T* update_u = house_mat + sub_length * sub_length;
274 T* mat_u_tmp = update_u + rows * rows;
276 for (
int i = 0; i < sub_length; ++i)
277 input_vec[i] = mat_b[(k + i) * cols + k];
280 &house_beta, epsilon, T(1));
282 house_beta, house_mat);
284 house_mat, sub_length, k, k);
286 for (
int i = k + 1; i < rows; ++i)
287 mat_b[i * cols + k] = T(0);
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)
295 for (
int j = 0; j < sub_length; ++j)
297 update_u[(k + i) * rows + (k + j)] =
298 house_mat[i * sub_length + j];
303 std::copy(mat_u, mat_u + rows * rows, mat_u_tmp);
310 for (
int i = k + 1; i < cols; ++i)
311 norm += mat_b[k * cols + i];
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
320 + inner_sub_length * inner_sub_length
321 + slice_rows * slice_cols
322 + slice_rows * slice_cols
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;
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;
335 for (
int i = 0; i < cols - k - 1; ++i)
336 inner_input_vec[i] = mat_b[k * cols + (k + 1 + i)];
339 inner_house_vec, &inner_house_beta, epsilon, norm);
341 inner_house_beta, inner_house_mat);
344 for (
int i = 0; i < slice_rows; ++i)
346 for (
int j = 0; j < slice_cols; ++j)
348 mat_b_tmp[i * slice_cols + j] =
349 mat_b[(k + i) * cols + (k + 1 + j)];
353 inner_house_mat, inner_sub_length, mat_b_res);
356 for (
int i = 0; i < slice_rows; ++i)
358 for (
int j = 0; j < slice_cols; ++j)
360 mat_b[(k + i) * cols + (k + 1 + j)] =
361 mat_b_res[i * slice_cols + j];
365 for (
int i = k + 2; i < cols; ++i)
366 mat_b[k * cols + i] = T(0);
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)
373 for (
int j = 0; j < inner_sub_length; ++j)
375 update_v[(k + i + 1) * cols + (k + j + 1)] =
376 inner_house_mat[i * inner_sub_length + j];
381 std::copy(mat_v, mat_v + cols * cols, mat_v_tmp);
393 int p,
int q, T
const& epsilon)
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;
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];
411 slice_length, mat_tmp);
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)];
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;
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)];
432 for (
int k = p; k < cols - q - 1; ++k)
434 T givens_c, givens_s;
441 alpha = mat_b[k * cols + k];
442 beta = mat_b[(k + 1) * cols + k];
444 internal::matrix_apply_givens_row(mat_b, cols, cols, k, k + 1,
446 internal::matrix_apply_givens_column(mat_q, rows, cols, k, k + 1,
449 if (k < (cols - q - 2))
451 alpha = mat_b[k * cols + (k + 1)];
452 beta = mat_b[k * cols + (k + 2)];
460 int row_index, T
const& epsilon)
462 for (
int i = row_index + 1; i < cols; ++i)
466 mat_b[row_index * cols + i] = T(0);
470 T norm =
MATH_POW2(mat_b[row_index * cols + i])
472 norm = std::sqrt(norm) *
MATH_SIGN(mat_b[i * cols + i]);
474 T givens_c = mat_b[i * cols + i] / norm;
475 T givens_s = mat_b[row_index * cols + i] / norm;
477 i, givens_c, givens_s);
479 i, givens_c, givens_s);
489 T* mat_u, T* vec_s, T* mat_v, T
const& epsilon)
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;
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;
507 mat_q_full, mat_b_full, mat_p, epsilon);
510 for (
int i = 0; i < rows; ++i)
512 for (
int j = 0; j < cols; ++j)
514 mat_q[i * cols + j] = mat_q_full[i * rows + j];
517 std::copy(mat_b_full, mat_b_full + cols * cols, mat_b);
520 int const max_iterations = rows * cols;
522 while (iteration < max_iterations)
527 for (
int i = 0; i < cols * cols; ++i)
529 T
const& entry = mat_b[i];
535 for (
int i = 0; i < (cols - 1); ++i)
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)]))
540 mat_b[i * cols + (i + 1)] = T(0);
547 for (
int k = 0; k < cols; ++k)
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)
553 for (
int j = 0; j < slice_len; ++j)
555 mat_b33[i * slice_len + j]
556 = mat_b[(cols - k - 1 + i) * cols + (cols - k - 1 + j)];
565 mat_b, cols, k + 1, epsilon))
579 T* mat_b22_tmp =
new T[(cols - q) * (cols - q)];
580 for (
int k = 0; k < (cols - q); ++k)
582 int const slice_len = k + 1;
583 for (
int i = 0; i < slice_len; ++i)
585 for (
int j = 0; j < slice_len; ++j)
587 mat_b22_tmp[i * slice_len + j]
588 = mat_b[(cols - q - k - 1 + i)
589 * cols + (cols - q - k - 1 + j)];
593 mat_b22_tmp, slice_len, slice_len, epsilon))
598 delete[] mat_b22_tmp;
600 int const p = cols - q - z;
606 bool diagonal_non_zero =
true;
608 for (nz = p; nz < (cols - q - 1); ++nz)
612 diagonal_non_zero =
false;
613 mat_b[nz * cols + nz] = T(0);
618 if (diagonal_non_zero)
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];
631 for (
int i = 0; i < cols; ++i)
633 if (vec_s[i] < epsilon)
635 vec_s[i] = -vec_s[i];
636 for (
int j = 0; j < rows; ++j)
638 int index = j * cols + i;
639 mat_u[index] = -mat_u[index];
652 T* mat_u, T* vec_s, T* mat_v, T
const& epsilon)
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;
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));
681 for (
int i = 0; i < length; ++i)
682 if (values[i] > largest)
700 T* mat_u, T* vec_s, T* mat_v, T
const& epsilon)
703 std::vector<T> mat_u_tmp;
704 std::vector<T> vec_s_tmp;
705 if (vec_s ==
nullptr)
707 vec_s_tmp.resize(cols);
708 vec_s = &vec_s_tmp[0];
710 std::vector<T> mat_v_tmp;
711 if (mat_v ==
nullptr)
713 mat_v_tmp.resize(cols * cols);
714 mat_v = &mat_v_tmp[0];
725 if (mat_u ==
nullptr)
727 mat_u_tmp.resize(rows * cols);
728 mat_u = &mat_u_tmp[0];
732 if (rows >= 5 * cols / 3)
734 internal::matrix_r_svd(mat_a, rows, cols,
735 mat_u, vec_s, mat_v, epsilon);
739 internal::matrix_gk_svd(mat_a, rows, cols,
740 mat_u, vec_s, mat_v, epsilon);
746 std::vector<T> mat_a_tmp(cols * cols, T(0));
747 std::copy(mat_a, mat_a + cols * rows, &mat_a_tmp[0]);
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);
755 if (mat_u !=
nullptr)
756 std::copy(&mat_u_tmp[0], &mat_u_tmp[0] + rows * cols, mat_u);
758 mat_u = &mat_u_tmp[0];
762 for (
int i = 0; i < cols; ++i)
764 int pos = internal::find_largest_ev_index(vec_s + i, cols - i);
778matrix_svd (T
const* mat_a,
int rows,
int cols,
779 T* mat_u, T* vec_s, T* mat_v, T
const& epsilon)
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)
787 mat_u_tmp.resize(rows * cols);
788 mat_u = &mat_u_tmp[0];
790 if (vec_s ==
nullptr)
792 vec_s_tmp.resize(cols);
793 vec_s = &vec_s_tmp[0];
795 if (mat_v ==
nullptr)
797 mat_v_tmp.resize(cols * cols);
798 mat_v = &mat_v_tmp[0];
811 std::vector<T> mat_a_tmp;
812 bool was_transposed =
false;
815 mat_a_tmp.resize(rows * cols);
816 std::copy(mat_a, mat_a + rows * cols, mat_a_tmp.begin());
820 mat_a = &mat_a_tmp[0];
821 was_transposed =
true;
825 if (rows >= 5 * cols / 3)
827 internal::matrix_r_svd(mat_a, rows, cols,
828 mat_u, vec_s, mat_v, epsilon);
832 internal::matrix_gk_svd(mat_a, rows, cols,
833 mat_u, vec_s, mat_v, epsilon);
842 for (
int i = rows; i < cols; ++i)
846 for (
int i = rows * rows - 1, y = rows - 1; y >= 0; --y)
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);
855 for (
int i = cols * rows - 1, y = cols - 1; y >= 0; --y)
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);
865 for (
int i = 0; i < cols; ++i)
867 int pos = internal::find_largest_ev_index(vec_s + i, cols - i);
879template <
typename T,
int M,
int N>
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;
888 matrix_svd<T>(mat_a.
begin(), M, N,
889 mat_u_ptr, mat_s_ptr, mat_v_ptr, epsilon);
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)
900template <
typename T,
int M,
int N>
911 for (
int i = 0; i < N; ++i)
915 S(i, i) = T(1) / S(i, i);
Matrix class for arbitrary dimensions and types.
Matrix< T, M, N > transposed(void) const
Returns a transposed copy of self by treating rows as columns.
#define MATH_INTERNAL_NAMESPACE_BEGIN
#define MATH_EPSILON_EQ(x, v, eps)
#define MATH_NAMESPACE_BEGIN
#define MATH_NAMESPACE_END
#define MATH_INTERNAL_NAMESPACE_END
#define MATH_SVD_DEFAULT_ZERO_THRESHOLD
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...
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.
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.
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.
void matrix_2x2_eigenvalues(T const *mat, T *smaller_ev, T *larger_ev)
Returns the larger eigenvalue of the given 2x2 matrix.
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,...
int find_largest_ev_index(T const *values, int length)
Returns the index of the largest eigenvalue.
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...
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 ...
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^...
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...
void matrix_svd_clear_super_entry(int rows, int cols, T *mat_b, T *mat_q, int row_index, T const &epsilon)
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...
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.
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...
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.
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.
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).
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.