feat: add SVD

tiankaima · Jan 3, 2024 · 52800ea · 52800ea
1 parent 0d4a132
commit 52800ea
Show file tree

Hide file tree

Showing 5 changed files with 162 additions and 7 deletions.
diff --git a/lib/IterationMethod/IterationMethod.h b/lib/IterationMethod/IterationMethod.h
@@ -7,8 +7,6 @@
 
 #include "../../includes/NLAMethods.h"
 
-#define ITERATION_METHOD_MAX_ITERATION 100000
-
 typedef struct {
     const Matrix &A;
     const Vector &b;

diff --git a/lib/SVDMethod/SVDMethod.cpp b/lib/SVDMethod/SVDMethod.cpp
@@ -4,6 +4,67 @@
 
 #include "SVDMethod.h"
 
+typedef struct {
+    Matrix R;
+    Matrix S;
+} PolarDecomposition2D_Result;
+
+PolarDecomposition2D_Result PolarDecomposition2D(const Matrix &A) {
+    auto x = A.matrix[0][0] + A.matrix[1][1];
+    auto y = A.matrix[1][0] - A.matrix[0][1];
+    auto r = std::sqrt(x * x + y * y);
+    auto c = x / r;
+    auto s = y / r;
+    auto R = Matrix(std::vector<std::vector<lld>>{
+            {c,  -s},
+            {s, c}
+    });
+    auto S = R.transpose() * A;
+    return {R, S};
+}
+
+WilkinsonShift_Result WilkinsonShiftIteration2D(const Matrix &A) {
+    auto [R, S] = PolarDecomposition2D(A);
+    auto V = Matrix::identity(2);
+    lld c, s, sigma_1, sigma_2;
+
+    if (S.matrix[0][1] == 0) {
+        c = 1;
+        s = 0;
+        sigma_1 = S.matrix[0][0];
+        sigma_2 = S.matrix[1][1];
+    } else {
+        auto tau = (S.matrix[0][0] - S.matrix[1][1]) / 2;
+        auto omega = std::sqrt(tau * tau + S.matrix[0][1] * S.matrix[0][1]);
+        auto t = S.matrix[0][1] / (tau + SIGN(tau) * omega);
+        c = 1 / std::sqrt(1 + t * t);
+        s = -c * t;
+        sigma_1 = S.matrix[0][0] * c * c + S.matrix[1][1] * s * s - 2 * S.matrix[0][1] * c * s;
+        sigma_2 = S.matrix[0][0] * s * s + S.matrix[1][1] * c * c + 2 * S.matrix[0][1] * c * s;
+    }
+
+    if (sigma_1 < sigma_2) {
+        std::swap(sigma_1, sigma_2);
+        V = Matrix(std::vector<std::vector<lld>>{
+                {-s, c},
+                {-c, -s}
+        });
+    } else {
+        V = Matrix(std::vector<std::vector<lld>>{
+                {c,  s},
+                {-s, c}
+        });
+    }
+
+    auto U = R * V;
+
+    auto B = Matrix(std::vector<std::vector<lld>>{
+            {sigma_1, 0},
+            {0,       sigma_2}
+    });
+
+    return {B, U.transpose(), V};
+}
 
 // SVD iteration with Wilkinson shift
 // W(A) -> P, Q, B, s.t. B = P * A * Q
@@ -14,6 +75,10 @@ WilkinsonShift_Result WilkinsonShiftIteration(const Matrix &matrix) {
     auto P = Matrix::identity(m);
     auto Q = Matrix::identity(n);
 
+    if (n == 2) {
+        return WilkinsonShiftIteration2D(matrix);
+    }
+
 #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);
@@ -79,6 +144,82 @@ ReformBidiagonalization_Result ReformBidiagonalization(const Matrix &matrix, ull
 }
 
 SVDMethod_Result SVDMethod(const Matrix &matrix) {
+#if DEBUG
+    if (matrix.rows < matrix.cols) {
+        throw std::invalid_argument("SVDMethod requires rows >= cols");
+    }
+#endif
+    auto [B, U, V] = BidiagonalizationMethod(matrix); // B = U * A * V
+    auto m = B.rows;
+    auto n = B.cols;
 
 
+    for (int count = 0; count < ITERATION_METHOD_MAX_ITERATION; count++) {
+        // clean B:
+        // TODO: Notice that these impl are not numerically stable, should be set to relative error
+        for (ull i = 0; i < n - 1; i++) {
+            if (std::abs(B.matrix[i][i + 1]) < 1e-8) {
+                B.matrix[i][i + 1] = 0;
+            }
+        }
+
+        for (ull i = 0; i < n; i++) {
+            if (std::abs(B.matrix[i][i]) < 1e-8) {
+                B.matrix[i][i] = 0;
+            }
+        }
+
+        ull pivot_i = 0;
+        ull pivot_j = 0;
+
+        for (ull i = n - 1; i > 0; i--) {
+            if (B.matrix[i - 1][i] != 0) {
+                pivot_j = i + 1;
+                break;
+            }
+        }
+
+        if (pivot_j == 0) {
+            return {B, U, V};
+        }
+
+        for (ull i = pivot_j - 1; i > 0; i--) {
+            if (B.matrix[i - 1][i] == 0) {
+                pivot_i = i;
+                break;
+            }
+        }
+
+        auto B22 = B.sub_matrix(pivot_i, pivot_j, pivot_i, pivot_j);
+
+        ull k = 0;
+        bool flag = false;
+        // iterate over B22, if there's a B22_ii = 0, use ReformBidiagonalization:
+        for (ull i = 0; i < B22.rows; i++) {
+            if (B22.matrix[i][i] == 0) {
+                flag = true;
+                k = i;
+                break;
+            }
+        }
+
+        if (flag) {
+            auto r = ReformBidiagonalization(B22, k);
+            B.set(pivot_i, pivot_j, pivot_j, pivot_j, r.B);
+            auto G_expanded = Matrix::identity(m);
+            G_expanded.set(pivot_i, pivot_j, pivot_i, pivot_j, r.G);
+            U = G_expanded * U ;
+        } else {
+            auto r = WilkinsonShiftIteration(B22);
+            B.set(pivot_i, pivot_j, pivot_i, pivot_j, r.B);
+
+            auto U_expanded = Matrix::identity(m);
+            U_expanded.set(pivot_i, pivot_j, pivot_i, pivot_j, r.P);
+            U = U_expanded * U;
+
+            auto V_expanded = Matrix::identity(n);
+            V_expanded.set(pivot_i, pivot_j, pivot_i, pivot_j, r.Q);
+            V = V * V_expanded;
+        }
+    }
 }
diff --git a/lib/SVDMethod/SVDMethod.h b/lib/SVDMethod/SVDMethod.h
@@ -8,18 +8,20 @@
 #include "../../includes/NLAMethods.h"
 
 typedef struct {
+    Matrix B;
     Matrix U;
-    Matrix S;
     Matrix V;
 } SVDMethod_Result;
 
+SVDMethod_Result SVDMethod(const Matrix &matrix);
 
 typedef struct {
     Matrix B;
     Matrix P;
     Matrix Q;
 } WilkinsonShift_Result;
 
+WilkinsonShift_Result WilkinsonShiftIteration2D(const Matrix &A);
 WilkinsonShift_Result WilkinsonShiftIteration(const Matrix &matrix);
 
 typedef struct {

diff --git a/lib/base.h b/lib/base.h
@@ -6,7 +6,7 @@
 #define NUMERICAL_ALGEBRA_BASE_H
 
 #define ENABLE_TIMING
-#define ITERATION_METHOD_MAX_ITERATION 100000
+#define ITERATION_METHOD_MAX_ITERATION 10000000
 
 #include "iostream"
 #include "iomanip"

diff --git a/main.cpp b/main.cpp
@@ -2,12 +2,26 @@
 
 // homework 8 (SVD):
 int main() {
-    auto A = Matrix(8, 8, 1, 10);
-    auto r = BidiagonalizationMethod(A);
+    auto A = Matrix(10, 8, 1, 10);
+    auto r = SVDMethod(A);
+
+    A.print();
 
     r.B.print();
 
-    (r.U * A * r.V).print();
+    (r.U * A * r.V).clean(1e-5).print();
+
+    (r.U.transpose() * r.U).print();
+    (r.V.transpose() * r.V).print();
 
     return 0;
+
+//    auto A = Matrix("[1 2; 0 4]");
+//
+//    auto r = WilkinsonShiftIteration2D(A);
+//
+//    r.B.print();
+//    (r.P.transpose() * A * r.Q.transpose()).print();
+//
+//    return 0;
 }