| 56 | |
| 57 | template <class T> |
| 58 | void SVD(const size_t dim[2], T* U_, T* S_, T* V_, T eps = -1) { |
| 59 | megdnn_assert(dim[0] >= dim[1]); |
| 60 | |
| 61 | { // Bi-diagonalization |
| 62 | size_t n = std::min(dim[0], dim[1]); |
| 63 | std::vector<T> house_vec(std::max(dim[0], dim[1])); |
| 64 | for (size_t i = 0; i < n; i++) { |
| 65 | // Column Householder |
| 66 | { |
| 67 | T x1 = S(i, i); |
| 68 | if (x1 < 0) |
| 69 | x1 = -x1; |
| 70 | |
| 71 | T x_inv_norm = 0; |
| 72 | for (size_t j = i; j < dim[0]; j++) { |
| 73 | x_inv_norm += S(j, i) * S(j, i); |
| 74 | } |
| 75 | if (x_inv_norm > 1e-7) |
| 76 | x_inv_norm = 1 / sqrt(x_inv_norm); |
| 77 | |
| 78 | T alpha = sqrt(1 + x1 * x_inv_norm); |
| 79 | T beta = x_inv_norm / alpha; |
| 80 | if (fabs(x_inv_norm) < 1e-7) |
| 81 | alpha = 0; // nothing to do |
| 82 | |
| 83 | house_vec[i] = -alpha; |
| 84 | for (size_t j = i + 1; j < dim[0]; j++) { |
| 85 | house_vec[j] = -beta * S(j, i); |
| 86 | } |
| 87 | if (S(i, i) < 0) |
| 88 | for (size_t j = i + 1; j < dim[0]; j++) { |
| 89 | house_vec[j] = -house_vec[j]; |
| 90 | } |
| 91 | } |
| 92 | for (size_t k = i; k < dim[1]; k++) { |
| 93 | T dot_prod = 0; |
| 94 | for (size_t j = i; j < dim[0]; j++) { |
| 95 | dot_prod += S(j, k) * house_vec[j]; |
| 96 | } |
| 97 | for (size_t j = i; j < dim[0]; j++) { |
| 98 | S(j, k) -= dot_prod * house_vec[j]; |
| 99 | } |
| 100 | } |
| 101 | for (size_t k = 0; k < dim[0]; k++) { |
| 102 | T dot_prod = 0; |
| 103 | for (size_t j = i; j < dim[0]; j++) { |
| 104 | dot_prod += U(k, j) * house_vec[j]; |
| 105 | } |
| 106 | for (size_t j = i; j < dim[0]; j++) { |
| 107 | U(k, j) -= dot_prod * house_vec[j]; |
| 108 | } |
| 109 | } |
| 110 | |
| 111 | // Row Householder |
| 112 | if (i >= n - 1) |
| 113 | continue; |
| 114 | { |
| 115 | T x1 = S(i, i + 1); |