// SVD.cpp : Defines the entry point for the console application. // #include "stdafx.h" // JordanGauss.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include #include #include #include #include #include #include using namespace boost::multiprecision; typedef number, et_off> cpp_bin_float_4096; using namespace std; std::vector > matrix_i; template void print_matrix(std::vector > matrix) { for (std::size_t row = 0; row < matrix.size(); row++) { for (std::size_t col = 0; col < matrix[row].size(); col++) std::cout << std::setprecision(1) << std::fixed << matrix[row][col] << " "; std::cout << "\n"; } //std::cout << "\n"; } template void print_vector(std::vector a_vector) { //############################################################## for (auto i = a_vector.begin(); i != a_vector.end(); ++i) std::cout << *i << ' '; std::cout << "\n"; } template void matrix_transpose(std::vector > matrix1, std::vector >& matrix2) { matrix2.resize(matrix1.size()); for (std::size_t row = 0; row < matrix1.size(); row++) { matrix2[row].resize(matrix1[row].size()); for (std::size_t col = 0; col < matrix1[row].size(); col++) matrix2[row][col] = matrix1[col][row]; } //############################################################## std::cout << "\n[matrix]transposed = \n"; print_matrix(matrix2); //############################################################## } template void matrix_by_matrix(std::vector > matrix1, std::vector >& matrix2, std::vector >& matrix3) { matrix3.resize(matrix1.size()); for (std::size_t row = 0; row < matrix1.size(); row++) { matrix3[row].resize(matrix1[row].size()); for (std::size_t col = 0; col < matrix1[row].size(); col++) { matrix3[row][col] = 0.00; for (std::size_t k = 0; k < matrix1[row].size(); k++) matrix3[row][col] += matrix1[row][k] * matrix2[k][col]; } } //############################################################## std::cout << "\n[matrices]multiplied = \n"; print_matrix(matrix3); //############################################################## } template void get_hermitian_matrix(std::vector eigenvector, std::vector >& h_matrix) { h_matrix.resize(eigenvector.size()); for (std::size_t row = 0; row < eigenvector.size(); row++) h_matrix[row].resize(eigenvector.size()); //----------------------------------------------------- // added to zero-out the matrix to avoid nan for (std::size_t row = 0; row < h_matrix.size(); row++) { for (std::size_t col = 0; col < h_matrix[row].size(); col++) h_matrix[row][col] = 0; } //----------------------------------------------------- h_matrix[0][0] = 1.0 / eigenvector[0]; for (std::size_t row = 1; row < eigenvector.size(); row++) h_matrix[row][0] = -eigenvector[row] / eigenvector[0]; for (std::size_t row = 1; row < eigenvector.size(); row++) h_matrix[row][row] = 1; //############################################################## std::cout << "\nHermitian matrix = \n"; print_matrix(h_matrix); //############################################################## } template void get_hermitian_matrix_inverse(std::vector eigenvector, std::vector >& ih_matrix) { ih_matrix.resize(eigenvector.size()); for (std::size_t row = 0; row < eigenvector.size(); row++) ih_matrix[row].resize(eigenvector.size()); //----------------------------------------------------- // added to zero-out the matrix to avoid nan for (std::size_t row = 0; row < ih_matrix.size(); row++) { for (std::size_t col = 0; col < ih_matrix[row].size(); col++) ih_matrix[row][col] = 0; } //----------------------------------------------------- ih_matrix[0][0] = eigenvector[0]; for (std::size_t row = 1; row < eigenvector.size(); row++) ih_matrix[row][0] = -eigenvector[row]; for (std::size_t row = 1; row < eigenvector.size(); row++) ih_matrix[row][row] = 1; //############################################################## std::cout << "\nHermitian matrix inverse = \n"; print_matrix(ih_matrix); //############################################################## } template void jordan_gaussian_transform( std::vector > matrix, std::vector& eigenvector) { const Arg eps = 10e-6; bool eigenv_found = false; for (std::size_t s = 0; s < matrix.size() - 1 && !eigenv_found; s++) { std::size_t col = s; Arg alpha = matrix[s][s]; while (col < matrix[s].size() && alpha != 0 && alpha != 1) matrix[s][col++] /= alpha; for (std::size_t col = s; col < matrix[s].size() && !alpha; col++) std::swap(matrix[s][col], matrix[s + 1][col]); for (std::size_t row = 0; row < matrix.size(); row++) { Arg gamma = matrix[row][s]; for (std::size_t col = s; col < matrix[row].size() && row != s; col++) matrix[row][col] = matrix[row][col] - matrix[s][col] * gamma; } } std::size_t row = 0; while (row < matrix.size()) eigenvector.push_back(-matrix[row++][matrix.size() - 1]); eigenvector[eigenvector.size() - 1] = 1; //############################################################## std::cout << "\nJordan-Gaussian transform: \n"; std::cout << "\nEigenvalues = \n"; print_vector(eigenvector); //############################################################## } template void get_inverse_diagonal_matrix(std::vector > matrix, std::vector >& inv_matrix) { inv_matrix.resize(matrix.size()); //----------------------------------------------------- // added to zero-out the matrix to avoid nan for (std::size_t row = 0; row < inv_matrix.size(); row++) { for (std::size_t col = 0; col < inv_matrix[row].size(); col++) inv_matrix[row][col] = 0; } //----------------------------------------------------- for (std::size_t index = 0; index < matrix.size(); index++) { inv_matrix[index].resize(matrix[index].size(), 0); inv_matrix[index][index] = 1.0 / matrix[index][index]; } //############################################################## std::cout << "\ninverse diagonal matrix = \n"; print_matrix(inv_matrix); //############################################################## } template void get_reduced_matrix(std::vector > matrix, std::vector >& r_matrix, std::size_t new_size) { if (new_size > 1) { r_matrix.resize(new_size); std::size_t index_d = matrix.size() - new_size; std::size_t row = index_d, row_n = 0; while (row < matrix.size()) { r_matrix[row_n].resize(new_size); std::size_t col = index_d, col_n = 0; while (col < matrix.size()) r_matrix[row_n][col_n++] = matrix[row][col++]; row++; row_n++; } } else if (new_size == 1) { r_matrix.resize(new_size); r_matrix[0].resize(new_size); r_matrix[0][0] = matrix[1][1]; } //############################################################## std::cout << "\nReduced matrix = \n"; print_matrix(r_matrix); //############################################################## } template void generate_matrix(std::vector >& \ matrix, std::size_t rows, std::size_t cols) { std::srand((unsigned int)std::time(nullptr)); matrix.resize(rows); for (std::size_t row = 0; row < matrix.size(); row++) { matrix[row].resize(cols); for (std::size_t col = 0; col < matrix[row].size(); col++) matrix[row][col] = std::rand() % 10; } } template void compute_evd(std::vector > matrix, std::vector& eigenvalues, std::vector >& eigenvectors, std::size_t eig_count) { std::size_t m_size = matrix.size(); std::vector vec; vec.resize(m_size); std::generate(vec.begin(), vec.end(), []() { return rand() % 100; }); if (eigenvalues.size() == 0 && eigenvectors.size() == 0) { eigenvalues.resize(m_size); eigenvectors.resize(eigenvalues.size()); matrix_i = matrix; } std::vector > m; m.resize(m_size); for (std::size_t row = 0; row < m_size; row++) m[row].resize(100); Arg lambda_old = 0; std::size_t index = 0; bool is_eval = false; while (is_eval == false) { for (std::size_t row = 0; row < m_size && (index % 100) == 0; row++) m[row].resize(m[row].size() + 100); for (std::size_t row = 0; row < m_size; row++) { m[row][index] = 0; for (std::size_t col = 0; col < m_size; col++) m[row][index] += matrix[row][col] * vec[col]; } for (std::size_t col = 0; col < m_size; col++) vec[col] = m[col][index]; if (index > 0) { Arg lambda = ((index - 1) > 0) ? \ (m[0][index] / m[0][index - 1]) : m[0][index]; Arg Change = boost::multiprecision::fabs(lambda - lambda_old); is_eval = (Change < /*10e-15*/10e-10) ? true : false; std::cout << "\nFor Index = "; std::cout << std::setprecision(0) << std::fixed << index << " "; std::cout << " ... "; std::cout << "Change in lambda = "; std::cout << std::setprecision(11) << std::fixed << Change << " "; std::cout << "\n"; eigenvalues[eig_count] = lambda; lambda_old = lambda; } index++; } std::vector > matrix_new; if (m_size > 1) { std::vector > matrix_target; matrix_target.resize(m_size); for (std::size_t row = 0; row < m_size; row++) { matrix_target[row].resize(m_size); for (std::size_t col = 0; col < m_size; col++) matrix_target[row][col] = (row == col) ? \ (matrix[row][col] - eigenvalues[eig_count]) : matrix[row][col]; } std::vector eigenvector; jordan_gaussian_transform(matrix_target, eigenvector); std::vector > hermitian_matrix; get_hermitian_matrix(eigenvector, hermitian_matrix); std::vector > ha_matrix_product; matrix_by_matrix(hermitian_matrix, matrix, ha_matrix_product); std::vector > inverse_hermitian_matrix; get_hermitian_matrix_inverse(eigenvector, inverse_hermitian_matrix); std::vector > iha_matrix_product; matrix_by_matrix(ha_matrix_product, inverse_hermitian_matrix, iha_matrix_product); get_reduced_matrix(iha_matrix_product, matrix_new, m_size - 1); } if (m_size <= 1) { for (std::size_t index = 0; index < eigenvalues.size(); index++) { Arg lambda = eigenvalues[index]; std::vector > matrix_target; matrix_target.resize(matrix_i.size()); for (std::size_t row = 0; row < matrix_i.size(); row++) { matrix_target[row].resize(matrix_i.size()); for (std::size_t col = 0; col < matrix_i.size(); col++) matrix_target[row][col] = (row == col) ? \ (matrix_i[row][col] - lambda) : matrix_i[row][col]; } eigenvectors.resize(matrix_i.size()); jordan_gaussian_transform(matrix_target, eigenvectors[index]); Arg eigsum_sq = 0; for (std::size_t v = 0; v < eigenvectors[index].size(); v++) eigsum_sq += boost::multiprecision::pow(eigenvectors[index][v], 2.0); for (std::size_t v = 0; v < eigenvectors[index].size(); v++) eigenvectors[index][v] /= boost::multiprecision::sqrt(eigsum_sq); eigenvalues[index] = boost::multiprecision::sqrt(eigenvalues[index]); } return; } compute_evd(matrix_new, eigenvalues, eigenvectors, eig_count + 1); return; } template void svd(std::vector > matrix, std::vector >& s, std::vector >& u, std::vector >& v) { std::vector > matrix_t; matrix_transpose(matrix, matrix_t); std::vector > matrix_product1; matrix_by_matrix(matrix, matrix_t, matrix_product1); std::vector > matrix_product2; matrix_by_matrix(matrix_t, matrix, matrix_product2); std::vector > u_1; std::vector > v_1; std::vector eigenvalues; compute_evd(matrix_product2, eigenvalues, v_1, 0); matrix_transpose(v_1, v); s.resize(matrix.size()); for (std::size_t index = 0; index < eigenvalues.size(); index++) { s[index].resize(eigenvalues.size(),0); s[index][index] = eigenvalues[index]; } std::vector > s_inverse; get_inverse_diagonal_matrix(s, s_inverse); std::vector > av_matrix; matrix_by_matrix(matrix, v, av_matrix); matrix_by_matrix(av_matrix, s_inverse, u); } int main(int argc, char* argv[]) { std::size_t matrix_size = 0; std::vector > u, v; std::vector > matrix, s; std::cout << "Singular Value Decomposition (SVD):\n\n"; std::cout << "Enter size of matrix N = (50x50 max): "; std::cin >> matrix_size; if (matrix_size <= 50) { generate_matrix(matrix, matrix_size, matrix_size); std::cout << "\nA = \n"; print_matrix(matrix); svd(matrix, s, u, v); std::cout << "\nS = \n"; print_matrix(s); std::cout << "\nU = \n"; print_matrix(u); std::cout << "\nV = \n"; print_matrix(v); } else std::cout << "Wrong matrix size... (matrix decomposition recommended)"; std::cin.get(); //std::cin.get(); return 0; }