From 52800eabe027fe92edd8390aa3fc3dd623cbf49e Mon Sep 17 00:00:00 2001 From: TianKai Ma Date: Wed, 3 Jan 2024 19:30:55 +0800 Subject: [PATCH] feat: add SVD --- lib/IterationMethod/IterationMethod.h | 2 - lib/SVDMethod/SVDMethod.cpp | 141 ++++++++++++++++++++++++++ lib/SVDMethod/SVDMethod.h | 4 +- lib/base.h | 2 +- main.cpp | 20 +++- 5 files changed, 162 insertions(+), 7 deletions(-) diff --git a/lib/IterationMethod/IterationMethod.h b/lib/IterationMethod/IterationMethod.h index b98bacf..fe2e063 100644 --- 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 index 1bbd306..43119f0 100644 --- 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>{ + {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>{ + {-s, c}, + {-c, -s} + }); + } else { + V = Matrix(std::vector>{ + {c, s}, + {-s, c} + }); + } + + auto U = R * V; + + auto B = Matrix(std::vector>{ + {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 index 5cc019f..ee24a8b 100644 --- a/lib/SVDMethod/SVDMethod.h +++ b/lib/SVDMethod/SVDMethod.h @@ -8,11 +8,12 @@ #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; @@ -20,6 +21,7 @@ typedef struct { 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 index abdcf23..9c0b860 100644 --- 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 index 5aadfa4..62217f1 100644 --- 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; } \ No newline at end of file