fix: Bidiagonalization algo, sync progress

lib/HouseholderMethod/Bidiagonalization.cpp

@@ -6,44 +6,57 @@
 
 Bidiagonalization_Result BidiagonalizationMethod(const Matrix &matrix) {
 #ifdef DEBUG
-    if (matrix.rows < matrix.cols) {
-        throw std::invalid_argument("BidiagonalizationMethod: rows < cols");
+    // TODO: I'm limiting rows > cols + 1 for now, should have been rows >= cols
+    if (matrix.rows <= matrix.cols + 1) {
+        throw std::invalid_argument("BidiagonalizationMethod requires rows >= cols");
     }
 #endif
 
-    auto A = matrix;
-    auto m = A.rows;
-    auto n = A.cols;
+    auto B = Matrix(matrix);
+    auto m = B.rows;
+    auto n = B.cols;
     auto U = Matrix::identity(m);
     auto V = Matrix::identity(n);
 
     for (ull k = 0; k < n; k++) {
-        auto x = A.sub_array_col(k, m, k);
-        auto [v, beta] = HouseHolderMethod(x);
+        auto [v, beta] = HouseHolderMethod(B.sub_array_col(k, m, k));
 
-        auto P = Matrix::identity(m - k);
-        P = P - product(v, v) * beta;
+        auto P = Matrix::identity(m);
+        auto P_sub = Matrix::identity(m - k) - product(v, v) * beta;
+        P.set(k, m, k, m, P_sub);
+        U = P * U;
 
-        auto P_extended = Matrix::identity(m);
-        P_extended.set(k, m, k, m, P);
+//        B = P * B;
 
-        A = P_extended.transpose() * A;
-        U = U * P_extended;
+        // 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 x = A.sub_array_row(k, k + 1, n);
-            auto [v, beta] = HouseHolderMethod(x);
-
-            auto Q = Matrix::identity(n - k - 1);
-            Q = Q - product(v, v) * beta;
-
-            auto Q_extended = Matrix::identity(n);
-            Q_extended.set(k + 1, n, k + 1, n, Q);
-
-            A = A * Q_extended;
-            V = V * Q_extended;
+            auto [v, beta] = HouseHolderMethod(B.sub_array_row(k, k + 1, n));
+
+            auto Q = Matrix::identity(n);
+            auto Q_sub = Matrix::identity(n - k - 1) - product(v, v) * beta;
+            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));
         }
     }
 
-    return {A, U, V};
+    return {B, U, V};
 }

lib/JacobiMethod/JacobiMethod.h

@@ -18,6 +18,8 @@ enum JacobiMethodType {
     JACOBI_METHOD_TYPE_CLASSIC, JACOBI_METHOD_TYPE_LOOP, JACOBI_METHOD_TYPE_THRESHOLD
 };
 
+Matrix RotationMatrix(ull n, ull p, ull q, lld c, lld s);
+
 /// P^T A P = H, returns H, P
 /// A must be symmetric
 JIterationMethodOutput JacobiMethod(const Matrix &matrix, JacobiMethodType type = JACOBI_METHOD_TYPE_CLASSIC);

lib/SVDMethod/SVDMethod.cpp

@@ -17,40 +17,47 @@ WilkinsonShift_Result WilkinsonShiftIteration(const Matrix &matrix) {
     auto P = Matrix::identity(n);
     auto Q = Matrix::identity(n);
 
-    // delta_n ^ 2 + gamma_(n-1) ^ 2
-    auto alpha = B.matrix[n - 1][n - 1] * B.matrix[n - 1][n - 1] + B.matrix[n - 2][n - 1] * B.matrix[n - 2][n - 1];
-    // (delta_(n-1) ^ 2 + gamma_(n-1) ^ 2 - alpha)/ 2
-    auto delta = (B.matrix[n - 2][n - 2] * B.matrix[n - 2][n - 2] + B.matrix[n - 2][n - 1] * B.matrix[n - 2][n - 1] -
-                  alpha) / 2;
-    // delta_(n-1) * gamma_(n-1)
-    auto beta = B.matrix[n - 2][n - 2] * B.matrix[n - 2][n - 1];
-
-    auto mu = alpha - beta * beta / (delta + SIGN(delta) * sqrt(delta * delta + beta * beta));
-    auto y = B.matrix[0][0] * B.matrix[0][0] - mu;
-    auto z = B.matrix[0][0] * B.matrix[0][1];
-
-    for (ull k = 0; k < n - 1; k++) {
-        auto tmpVector = Vector(std::vector<lld>{y, z});
-        auto [b, beta] = HouseHolderMethod(tmpVector);
-        auto w = product(b, b) * beta;
-
-        auto q = MAX(1, k);
-        auto B_sub = B.sub_matrix(k, k + 2, q - 1, n);
-        B_sub = B_sub - w * B_sub;
-        B.set(k, k + 2, q - 1, n, B_sub);
-
-        auto r = MIN(k + 2, n);
-        B_sub = B.sub_matrix(0, r, k, k + 2);
-        B_sub = B_sub - B_sub * w;
-        B.set(0, r, k, k + 2, B_sub);
-
-        auto P_sub = P.sub_matrix(0, n, k, k + 2);
-        P_sub = P_sub - P_sub * w;
-        P.set(0, n, k, k + 2, P_sub);
-
-        y = B.matrix[k + 1][k];
-        if (k < n - 2) {
-            z = B.matrix[k + 2][k];
+#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);
+
+    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);
+
+        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);
+
+        t = -z / y; // 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;
         }
     }
+
+    return {B, P, Q};
 }

main.cpp

@@ -1,74 +1,17 @@
 #include "includes/NLAMethods.h"
 
-Matrix A_7_1(int n) {
-    auto A = Matrix(n, n);
-    for (int i = 0; i < n; i++) {
-        for (int j = 0; j < n; j++) {
-            if (i == j) {
-                A.matrix[i][j] = 4;
-            }
-            if (i == j + 1 || i == j - 1) {
-                A.matrix[i][j] = 1;
-            }
-        }
-    }
-    return A;
-}
-
-Matrix A_7_2(int n) {
-    auto A = Matrix(n, n);
-    for (ull i = 0; i < n; i++) {
-        for (ull j = 0; j < n; j++) {
-            if (i == j) {
-                A.matrix[i][j] = 2;
-            } else if (i == j + 1 || i == j - 1) {
-                A.matrix[i][j] = -1;
-            }
-        }
-    }
-
-    return A;
-}
-
-void part_1() {
-    for (auto n: {50, 60, 70, 80, 90, 100}) {
-        std::cout << "n = " << n << std::endl;
-
-        auto A = A_7_1(n);
-        auto r = JacobiMethod(A, JACOBI_METHOD_TYPE_THRESHOLD);
-
-        r.result.A.clean(1e-8).print();
-        r.result.P.print();
-
-        std::cout << "Iteration count: " << r.iteration_count << std::endl;
-        std::cout << "Time cost: " << r.time_cost.count() << std::endl;
-    }
-}
-
-void part_2() {
-    auto n = 100;
-    auto x = Vector(n, 2);
-    auto y = Vector(n - 1, -1);
-    auto A = A_7_2(n);
+// homework 8 (SVD):
+int main() {
+    auto A = Matrix(5, 3, 1, 10);
+//    A.print();
 
+    auto r = BidiagonalizationMethod(A);
 
-    auto I = Matrix::identity(n);
-    auto x_default = Vector(n, 1);
-    auto r = BisectMethod(x, y);
+    r.B.print();
+//    r.U.print();
+//    r.V.print();
 
-    for (auto &lambda: r.array) {
-        auto B = A - I * lambda;
-        auto input = PowerIterationInput{B, x_default, 1000,};
-        auto k = RevPowerIteration(input);
-        std::cout << std::fixed << std::setprecision(6);
-        std::cout << "lambda = " << lambda << "\t :";
-        k.result.print();
-    }
-}
+    (r.U * A * r.V).print();
 
-// homework 7:
-int main() {
-    part_1();
-    part_2();
     return 0;
 }