Skip to content
This repository has been archived by the owner on Jun 21, 2024. It is now read-only.

Commit

Permalink
feat: finish homework 8 (SVD)
Browse files Browse the repository at this point in the history 
last programming homework, finally
additional credit goes to cffjiang@UCLA, whose work could be found here:

https://www.math.ucla.edu/~cffjiang/research/svd/svd.pdf

which is absolutely amazing, fixing the last puzzle of 2x2 SVD for me.

Also thx to TA & all teachers for their incredibly great help.
  • Loading branch information
@tiankaima
tiankaima committed Jan 3, 2024
1 parent f4f1c41 commit 068a107
Show file tree
Hide file tree
Showing 4 changed files with 69 additions and 21 deletions.
17 changes: 13 additions & 4 deletions lib/SVDMethod/SVDMethod.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ ReformBidiagonalization_Result ReformBidiagonalization(const Matrix &matrix, ull
return {B, G};
}

SVDMethod_Result SVDMethod(const Matrix &matrix) {
SVDMethod_Output SVDMethod(const Matrix &matrix) {
#if DEBUG
if (matrix.rows < matrix.cols) {
throw std::invalid_argument("SVDMethod requires rows >= cols");
Expand All @@ -153,18 +153,19 @@ SVDMethod_Result SVDMethod(const Matrix &matrix) {
auto m = B.rows;
auto n = B.cols;

TIMING_START

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) {
if (std::abs(B.matrix[i][i + 1]) < 1e-4) {
B.matrix[i][i + 1] = 0;
}
}

for (ull i = 0; i < n; i++) {
if (std::abs(B.matrix[i][i]) < 1e-8) {
if (std::abs(B.matrix[i][i]) < 1e-4) {
B.matrix[i][i] = 0;
}
}
Expand All @@ -180,7 +181,15 @@ SVDMethod_Result SVDMethod(const Matrix &matrix) {
}

if (pivot_j == 0) {
return {B, U, V};
TIMING_END

return {
B,
U,
V,
count,
TIMING_RETURN_DURATION
};
}

for (ull i = pivot_j - 1; i > 0; i--) {
Expand Down
4 changes: 3 additions & 1 deletion lib/SVDMethod/SVDMethod.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,9 @@ typedef struct {
Matrix V;
} SVDMethod_Result;

SVDMethod_Result SVDMethod(const Matrix &matrix);
using SVDMethod_Output = IterationMethodOutput<SVDMethod_Result>;

SVDMethod_Output SVDMethod(const Matrix &matrix);

typedef struct {
Matrix B;
Expand Down
1 change: 1 addition & 0 deletions lib/base.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "random"
#include "chrono"
#include "optional"
#include "algorithm"

#define ull unsigned long long
#define lld long double
Expand Down
68 changes: 52 additions & 16 deletions main.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,27 +1,63 @@
#include "includes/NLAMethods.h"

auto A = Matrix(
std::vector<std::vector<lld>>{
{1.0000000000, 4.9176000000, 1.0000000000, 3.4720000000, 0.9980000000, 1.0000000000, 7.0000000000, 4.0000000000, 42.0000000000, 3.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 5.0208000000, 1.0000000000, 3.5310000000, 1.5000000000, 2.0000000000, 7.0000000000, 4.0000000000, 62.0000000000, 1.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 4.5429000000, 1.0000000000, 2.2750000000, 1.1750000000, 1.0000000000, 6.0000000000, 3.0000000000, 40.0000000000, 2.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 4.5573000000, 1.0000000000, 4.0500000000, 1.2320000000, 1.0000000000, 6.0000000000, 3.0000000000, 54.0000000000, 4.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 5.0597000000, 1.0000000000, 4.4550000000, 1.1210000000, 1.0000000000, 6.0000000000, 3.0000000000, 42.0000000000, 3.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 3.8910000000, 1.0000000000, 4.4550000000, 0.9880000000, 1.0000000000, 6.0000000000, 3.0000000000, 56.0000000000, 2.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 5.8980000000, 1.0000000000, 5.8500000000, 1.2400000000, 1.0000000000, 7.0000000000, 3.0000000000, 51.0000000000, 2.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 5.6039000000, 1.0000000000, 9.5200000000, 1.5010000000, 0.0000000000, 6.0000000000, 3.0000000000, 32.0000000000, 1.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 15.4202000000, 2.5000000000, 9.8000000000, 3.4200000000, 2.0000000000, 10.0000000000, 5.0000000000, 42.0000000000, 2.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 14.4598000000, 2.5000000000, 12.8000000000, 3.0000000000, 2.0000000000, 9.0000000000, 5.0000000000, 14.0000000000, 4.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 5.8282000000, 1.0000000000, 6.4350000000, 1.2250000000, 2.0000000000, 6.0000000000, 3.0000000000, 32.0000000000, 1.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 5.3003000000, 1.0000000000, 4.9883000000, 1.5520000000, 1.0000000000, 6.0000000000, 3.0000000000, 30.0000000000, 1.0000000000, 2.0000000000, 0.0000000000},
{1.0000000000, 6.2712000000, 1.0000000000, 5.5200000000, 0.9750000000, 1.0000000000, 5.0000000000, 2.0000000000, 30.0000000000, 1.0000000000, 2.0000000000, 0.0000000000},
{1.0000000000, 5.9592000000, 1.0000000000, 6.6660000000, 1.1210000000, 2.0000000000, 6.0000000000, 3.0000000000, 32.0000000000, 2.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 5.0500000000, 1.0000000000, 5.0000000000, 1.0200000000, 0.0000000000, 5.0000000000, 2.0000000000, 46.0000000000, 4.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 5.6039000000, 1.0000000000, 9.5200000000, 1.5010000000, 0.0000000000, 6.0000000000, 3.0000000000, 32.0000000000, 1.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 8.2464000000, 1.5000000000, 5.1500000000, 1.6640000000, 2.0000000000, 8.0000000000, 4.0000000000, 50.0000000000, 4.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 6.6969000000, 1.5000000000, 6.0920000000, 1.4880000000, 1.5000000000, 7.0000000000, 3.0000000000, 22.0000000000, 1.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 7.7841000000, 1.5000000000, 7.1020000000, 1.3760000000, 1.0000000000, 6.0000000000, 3.0000000000, 17.0000000000, 2.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 9.0384000000, 1.0000000000, 7.8000000000, 1.5000000000, 1.5000000000, 7.0000000000, 3.0000000000, 23.0000000000, 3.0000000000, 3.0000000000, 0.0000000000},
{1.0000000000, 5.9894000000, 1.0000000000, 5.5200000000, 1.2560000000, 2.0000000000, 6.0000000000, 3.0000000000, 40.0000000000, 4.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 7.5422000000, 1.5000000000, 4.0000000000, 1.6900000000, 1.0000000000, 6.0000000000, 3.0000000000, 22.0000000000, 1.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 8.7951000000, 1.5000000000, 9.8900000000, 1.8200000000, 2.0000000000, 8.0000000000, 4.0000000000, 50.0000000000, 1.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 6.0931000000, 1.5000000000, 6.7265000000, 1.6520000000, 1.0000000000, 6.0000000000, 3.0000000000, 44.0000000000, 4.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 8.3607000000, 1.5000000000, 9.1500000000, 1.7770000000, 2.0000000000, 8.0000000000, 4.0000000000, 48.0000000000, 1.0000000000, 1.0000000000, 1.0000000000},
{1.0000000000, 8.1400000000, 1.0000000000, 8.0000000000, 1.5040000000, 2.0000000000, 7.0000000000, 3.0000000000, 3.0000000000, 1.0000000000, 3.0000000000, 0.0000000000},
{1.0000000000, 9.1416000000, 1.5000000000, 7.3262000000, 1.8310000000, 1.5000000000, 8.0000000000, 4.0000000000, 31.0000000000, 4.0000000000, 1.0000000000, 0.0000000000},
{1.0000000000, 12.0000000000, 1.5000000000, 5.0000000000, 1.2000000000, 2.0000000000, 6.0000000000, 3.0000000000, 30.0000000000, 3.0000000000, 1.0000000000, 1.0000000000}
}
);

// homework 8 (SVD):
int main() {
auto A = Matrix(10, 8, 1, 10);
auto r = SVDMethod(A);
auto m = A.rows;
auto n = A.cols;
auto [r, count, timing] = SVDMethod(A);

A.print();
std::cout << "iterations: " << count << std::endl;
std::cout << "time spent: " << timing.count() << "μs" << std::endl;

r.B.print();
auto eigenvalues = Vector(r.B.cols);
// take diagonal elements of B as eigenvalues:
for (ull i = 0; i < r.B.cols; i++) {
eigenvalues.array[i] = r.B.matrix[i][i];
}

(r.U * A * r.V).clean(1e-5).print();
// sort eigenvalues in descending order:
std::sort(eigenvalues.array.begin(), eigenvalues.array.end(), std::less<>());

(r.U.transpose() * r.U).print();
(r.V.transpose() * r.V).print();
// print eigenvalues:
std::cout << "eigenvalues = ";
eigenvalues.print();

return 0;
std::cout << "ep = " << MatrixNorm_Infinity(r.U.transpose() * r.U - Matrix::identity(m)) << std::endl
<< "eq = " << MatrixNorm_Infinity(r.V.transpose() * r.V - Matrix::identity(n)) << std::endl
<< "et = " << MatrixNorm_Infinity(r.U * A * r.V - r.B) << std::endl;

// auto A = Matrix("[1 2; 0 4]");
//
// auto r = WilkinsonShiftIteration2D(A);
//
// r.B.print();
// (r.P.transpose() * A * r.Q.transpose()).print();
//
// return 0;
return 0;
}

0 comments on commit 068a107

Please sign in to comment.