fix: WilkinsonShiftIteration

Author: tiankaima

lib/HouseholderMethod/Bidiagonalization.cpp

@@ -26,14 +26,9 @@ Bidiagonalization_Result BidiagonalizationMethod(const Matrix &matrix) {
         P.set(k, m, k, m, P_sub);
         U = P * U;
 
-//        B = P * B;
-
-        // TODO: WTF WHY DOESN'T THIS WORK FUCK IT
         auto A_sub = B.sub_matrix(k, m, k, n);
         A_sub = A_sub - product(v, v) * beta * A_sub;
         B.set(k, m, k, n, A_sub);
-//
-//        B.set_col(k + 1, m, k, v.sub_vector(1, m - k));
 
         if (k < n - 2) {
             auto [v, beta] = HouseHolderMethod(B.sub_array_row(k, k + 1, n));
@@ -43,18 +38,9 @@ Bidiagonalization_Result BidiagonalizationMethod(const Matrix &matrix) {
             Q.set(k + 1, n, k + 1, n, Q_sub);
             V = Q * V;
 
-//            B = B * Q;
-
-            // TODO: WTF WHY DOESN'T THIS WORK FUCK IT
-            // NEVER MIND, the B.set_row was meant to store Q,
-            // but let's skip it as it's really confusing
-            // > what a wierd way to store stuff,
-            // > some FORTRAN programing perhaps?
             auto A_sub = B.sub_matrix(k, m, k + 1, n);
             A_sub = A_sub - A_sub * product(v, v) * beta;
             B.set(k, m, k + 1, n, A_sub);
-//
-//            B.set_row(k, k + 2, n, v.sub_vector(1, n - k - 1));
         }
     }
 

lib/SVDMethod/SVDMethod.cpp

@@ -6,75 +6,75 @@
 
 
 // SVD iteration with Wilkinson shift
+// W(A) -> P, Q, B, s.t. B = P * A * Q
 WilkinsonShift_Result WilkinsonShiftIteration(const Matrix &matrix) {
     auto B = Matrix(matrix);
-    auto n = B.rows;
-    auto P = Matrix::identity(n);
+    auto m = B.rows;
+    auto n = B.cols;
+    auto P = Matrix::identity(m);
     auto Q = Matrix::identity(n);
 
 #define DELTA(i) B.matrix[i - 1][i - 1]
 #define GAMMA(i) B.matrix[i - 1][i]
-
     auto alpha = DELTA(n) * DELTA(n) + GAMMA(n - 1) * GAMMA(n - 1);
     auto delta = (DELTA(n - 1) * DELTA(n - 1) + GAMMA(n - 2) * GAMMA(n - 2) - alpha) / 2;
     auto beta = DELTA(n - 1) * GAMMA(n - 1);
     auto mu = alpha - beta * beta / (delta + SIGN(delta) * std::sqrt(delta * delta + beta * beta));
 
     auto y = DELTA(1) * DELTA(1) - mu;
     auto z = DELTA(1) * GAMMA(1);
+#undef DELTA
+#undef GAMMA
 
-    for (ull k = 0; k < n; k++) {
-        auto t = -z / y; // tan(theta)
-        auto c = 1 / std::sqrt(1 + t * t); // cos(theta)
-        auto s = c * t; // sin(theta)
-
-        Q = Q * RotationMatrix(n, k + 1, k + 2, c, s);
+    auto t = -z / y;
+    auto c = 1 / std::sqrt(1 + t * t);
+    auto s = c * t;
 
-        y = c * DELTA(k + 1) - s * GAMMA(k + 1);
-        z = -s * DELTA(k + 2);
-        GAMMA(k + 1) = s * DELTA(k + 1) + c * GAMMA(k + 1);
-        DELTA(k + 2) = c * DELTA(k + 2);
+    auto G1 = RotationMatrix(n, 0, 1, c, s);
+    B = B * G1;
+    Q = Q * G1;
 
-        t = -z / y; // tan(theta)
+    for (ull k = 0; k < n - 1; k++) {
+        t = -B.matrix[k + 1][k] / B.matrix[k][k]; // tan(theta)
         c = 1 / std::sqrt(1 + t * t); // cos(theta)
         s = c * t; // sin(theta)
 
-        P = RotationMatrix(n, k + 1, k + 2, c, s) * P;
-        if (k != n - 1) {
-            y = c * GAMMA(k + 1) - s * DELTA(k + 2);
-            z = -s * GAMMA(k + 2);
-            DELTA(k + 1) = s * GAMMA(k + 1) + c * DELTA(k + 2);
-            GAMMA(k + 2) = c * GAMMA(k + 2);
-        } else {
-            auto gamma_k = GAMMA(k + 1);
-            auto delta_k = DELTA(k + 2);
-
-            GAMMA(k + 1) = c * gamma_k - s * delta_k;
-            DELTA(k + 2) = s * gamma_k + c * delta_k;
+        auto G = RotationMatrix(m, k, k + 1, c, s);
+        B = G * B;
+        P = G * P;
+
+
+        if (k != n - 2) {
+            t = B.matrix[k][k + 2] / B.matrix[k][k + 1]; // tan(theta)
+            c = 1 / std::sqrt(1 + t * t); // cos(theta)
+            s = c * t; // sin(theta)
+
+            G = RotationMatrix(n, k + 1, k + 2, c, s);
+            B = B * G;
+            Q = Q * G;
         }
     }
 
     return {B, P, Q};
-#undef DELTA
-#undef GAMMA
 }
 
 // Consider position k have DELTA(k + 1) = 0
 // Then we can use Givens rotation to make GAMMA(k) = 0 as well:
 ReformBidiagonalization_Result ReformBidiagonalization(const Matrix &matrix, ull k) {
     auto B = Matrix(matrix);
-    auto n = B.rows;
-    auto G = Matrix::identity(n);
+    auto m = B.rows;
+    auto n = B.cols;
+    auto G = Matrix::identity(m);
 
     B.print();
 
-    for (ull i = k + 1; i < B.cols; i++) {
-        auto t = B.matrix[k][i] / B.matrix[i][i] ; // tan(theta)
+    for (ull i = k + 1; i < n; i++) {
+        auto t = B.matrix[k][i] / B.matrix[i][i]; // tan(theta)
         auto c = 1 / std::sqrt(1 + t * t); // cos(theta)
         auto s = c * t; // sin(theta)
 
-        G = RotationMatrix(n, k, i, c, s) * G;
-        B = RotationMatrix(n, k, i, c, s) * B;
+        G = RotationMatrix(m, k, i, c, s) * G;
+        B = RotationMatrix(m, k, i, c, s) * B;
 
         B.print();
         G.print();
@@ -83,6 +83,7 @@ ReformBidiagonalization_Result ReformBidiagonalization(const Matrix &matrix, ull
     return {B, G};
 }
 
-//SVDMethod_Result SVDMethod(const Matrix &matrix) {
-//
-//}
+SVDMethod_Result SVDMethod(const Matrix &matrix) {
+
+
+}

main.cpp

@@ -4,12 +4,14 @@
 int main() {
     auto A = Matrix(8, 5, 1, 10);
     A = BidiagonalizationMethod(A).B;
-    A.matrix[1][1] = 0;
+    A.matrix[1][0] = 5;
 
-    auto r = ReformBidiagonalization(A, 1);
+    A.print();
+
+    auto r = WilkinsonShiftIteration(A);
 
     r.B.print();
-    r.G.print();
+    (r.P * A * r.Q).print();
 
     return 0;
 }