Skip to content

Commit cca58e3

Browse files
kaydotdevsiriak
kaydotdev
and
siriak
authored
Add eigenvalue and eigenvector power iteration method (#177)
* Update vector and matrix extensions * Introduce fluent assertions Co-authored-by: Andrii Siriak <siryaka@gmail.com>
1 parent ca3e4bf commit cca58e3

File tree

15 files changed

+784
-503
lines changed

15 files changed

+784
-503
lines changed

‎Algorithms.Tests/Algorithms.Tests.csproj

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
<IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>
2020
<PrivateAssets>all</PrivateAssets>
2121
</PackageReference>
22+
<PackageReference Include="FluentAssertions" Version="5.10.3" />
2223
<PackageReference Include="Microsoft.NET.Test.Sdk" Version="16.3.0" />
2324
<PackageReference Include="nunit" Version="3.12.0" />
2425
<PackageReference Include="NUnit3TestAdapter" Version="3.15.1" />

‎Algorithms.Tests/LinearAlgebra/Eigenvalue/PowerIterationTests.cs

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
using System;
2+
using Algorithms.LinearAlgebra.Eigenvalue;
3+
using FluentAssertions;
4+
using Utilities.Extensions;
5+
using NUnit.Framework;
6+
7+
namespace Algorithms.Tests.LinearAlgebra.Eigenvalue
8+
{
9+
public class PowerIterationTests
10+
{
11+
private readonly double epsilon = Math.Pow(10, -5);
12+
13+
[Test]
14+
public void Dominant_ShouldThrowArgumentException_WhenSourceMatrixIsNotSquareShaped()
15+
{
16+
// Arrange
17+
var source = new double[,] {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}};
18+
19+
// Act
20+
Action action = () => PowerIteration.Dominant(source, StartVector(source.GetLength(0)), epsilon);
21+
22+
// Assert
23+
action.Should().Throw<ArgumentException>().WithMessage("The source matrix is not square-shaped.");
24+
}
25+
26+
[Test]
27+
public void Dominant_ShouldThrowArgumentException_WhenStartVectorIsNotSameSizeAsMatrix()
28+
{
29+
// Arrange
30+
var source = new double[,] {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
31+
var startVector = new double[] {1, 0, 0, 0};
32+
33+
// Act
34+
Action action = () => PowerIteration.Dominant(source, startVector, epsilon);
35+
36+
// Assert
37+
action.Should().Throw<ArgumentException>()
38+
.WithMessage("The length of the start vector doesn't equal the size of the source matrix.");
39+
}
40+
41+
[Test, TestCaseSource(nameof(DominantVectorTestCases))]
42+
public void Dominant_ShouldCalculateDominantEigenvalueAndEigenvector(
43+
double eigenvalue, double[] eigenvector, double[,] source)
44+
{
45+
// Act
46+
var (actualEigVal, actualEigVec) = PowerIteration.Dominant(source, StartVector(source.GetLength(0)), epsilon);
47+
48+
// Assert
49+
actualEigVal.Should().BeApproximately(eigenvalue, epsilon);
50+
actualEigVec.Magnitude().Should().BeApproximately(eigenvector.Magnitude(), epsilon);
51+
}
52+
53+
static readonly object[] DominantVectorTestCases =
54+
{
55+
new object[]
56+
{
57+
3.0,
58+
new[] { 0.7071039, 0.70710966 },
59+
new[,] { { 2.0, 1.0 }, { 1.0, 2.0 } }
60+
},
61+
new object[]
62+
{
63+
4.235889,
64+
new[] { 0.91287093, 0.40824829 },
65+
new[,] { { 2.0, 5.0 }, { 1.0, 2.0 } }
66+
}
67+
};
68+
69+
private double[] StartVector(int length) => new Random(Seed: 111111).NextVector(length);
70+
}
71+
}

‎Algorithms.Tests/Numeric/Decomposition/SVDTests.cs

Lines changed: 134 additions & 135 deletions
Original file line numberDiff line numberDiff line change
@@ -1,136 +1,135 @@
1-
using Algorithms.Numeric.Decomposition;
2-
using NUnit.Framework;
3-
using System;
4-
using M = Utilities.Extensions.MatrixExtensions;
5-
using V = Utilities.Extensions.VectorExtensions;
6-
7-
namespace Algorithms.Tests.Numeric.Decomposition
8-
{
9-
public class SvdTests
10-
{
11-
public void AssertMatrixEqual(double[,] matrix1, double[,] matrix2, double epsilon)
12-
{
13-
Assert.AreEqual(matrix1.GetLength(0), matrix2.GetLength(0));
14-
Assert.AreEqual(matrix1.GetLength(1), matrix2.GetLength(1));
15-
for (int i = 0; i < matrix1.GetLength(0); i++)
16-
{
17-
for (int j = 0; j < matrix1.GetLength(1); j++)
18-
{
19-
Assert.AreEqual(matrix1[i, j], matrix2[i, j], epsilon, $"At index ({i}, {j})");
20-
}
21-
}
22-
}
23-
24-
public double[,] GenerateRandomMatrix(int m, int n)
25-
{
26-
double[,] result = new double[m, n];
27-
Random random = new Random();
28-
for (int i = 0; i < m; i++)
29-
{
30-
for (int j = 0; j < n; j++)
31-
{
32-
result[i, j] = random.NextDouble() - 0.5;
33-
}
34-
}
35-
return result;
36-
}
37-
38-
[Test]
39-
public void RandomUnitVector()
40-
{
41-
double epsilon = 0.0001;
42-
// unit vector should have length 1
43-
Assert.AreEqual(1, V.Magnitude(ThinSvd.RandomUnitVector(10)), epsilon);
44-
// unit vector with single element should be [-1] or [+1]
45-
Assert.AreEqual(1, Math.Abs(ThinSvd.RandomUnitVector(1)[0]), epsilon);
46-
// two randomly generated unit vectors should not be equal
47-
Assert.AreNotEqual(ThinSvd.RandomUnitVector(10), ThinSvd.RandomUnitVector(10));
48-
}
49-
50-
[Test]
51-
public void Svd_Decompose()
52-
{
53-
CheckSvd(new double[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } });
54-
CheckSvd(new double[,] { { 1, 2, 3 }, { 4, 5, 6 } });
55-
CheckSvd(new double[,] { { 1, 0, 0, 0, 2 }, { 0, 3, 0, 0, 0 }, { 0, 0, 0, 0, 0 }, { 0, 2, 0, 0, 0 } });
56-
}
57-
58-
[Test]
59-
public void Svd_Random([Random(3, 10, 5)] int m, [Random(3, 10, 5)] int n)
60-
{
61-
double[,] matrix = GenerateRandomMatrix(m, n);
62-
CheckSvd(matrix);
63-
}
64-
65-
public void CheckSvd(double[,] testMatrix)
66-
{
67-
double epsilon = 1E-5;
68-
double[,] u;
69-
double[,] v;
70-
double[] s;
71-
(u, s, v) = ThinSvd.Decompose(testMatrix, 1E-8, 1000);
72-
73-
for (int i = 1; i < s.Length; i++)
74-
{
75-
// singular values should be arranged from greatest to smallest
76-
Assert.GreaterOrEqual(s[i - 1], s[i]);
77-
}
78-
79-
for (int i = 0; i < u.GetLength(1); i++)
80-
{
81-
double[] extracted = new double[u.GetLength(0)];
82-
// extract a column of u
83-
for (int j = 0; j < extracted.Length; j++)
84-
{
85-
extracted[j] = u[j, i];
86-
}
87-
88-
if (s[i] > epsilon)
89-
{
90-
// if the singular value is non-zero, then the basis vector in u should be a unit vector
91-
Assert.AreEqual(1, V.Magnitude(extracted), epsilon);
92-
}
93-
else
94-
{
95-
// if the singular value is zero, then the basis vector in u should be zeroed out
96-
Assert.AreEqual(0, V.Magnitude(extracted), epsilon);
97-
}
98-
}
99-
100-
for (int i = 0; i < v.GetLength(1); i++)
101-
{
102-
double[] extracted = new double[v.GetLength(0)];
103-
// extract column of v
104-
for (int j = 0; j < extracted.Length; j++)
105-
{
106-
extracted[j] = v[j, i];
107-
}
108-
109-
if (s[i] > epsilon)
110-
{
111-
// if the singular value is non-zero, then the basis vector in v should be a unit vector
112-
Assert.AreEqual(1, V.Magnitude(extracted), epsilon);
113-
}
114-
else
115-
{
116-
// if the singular value is zero, then the basis vector in v should be zeroed out
117-
Assert.AreEqual(0, V.Magnitude(extracted), epsilon);
118-
}
119-
}
120-
121-
// convert singular values to a diagonal matrix
122-
double[,] expanded = new double[s.Length, s.Length];
123-
for (int i = 0; i < s.Length; i++)
124-
{
125-
expanded[i, i] = s[i];
126-
}
127-
128-
129-
// matrix = U * S * V^t, definition of Singular Vector Decomposition
130-
AssertMatrixEqual(testMatrix,
131-
M.MultiplyGeneral(M.MultiplyGeneral(u, expanded), M.Transpose(v)), epsilon);
132-
AssertMatrixEqual(testMatrix,
133-
M.MultiplyGeneral(u, M.MultiplyGeneral(expanded, M.Transpose(v))), epsilon);
134-
}
135-
}
1+
using Algorithms.Numeric.Decomposition;
2+
using NUnit.Framework;
3+
using System;
4+
using Utilities.Extensions;
5+
using M = Utilities.Extensions.MatrixExtensions;
6+
using V = Utilities.Extensions.VectorExtensions;
7+
8+
namespace Algorithms.Tests.Numeric.Decomposition
9+
{
10+
public class SvdTests
11+
{
12+
public void AssertMatrixEqual(double[,] matrix1, double[,] matrix2, double epsilon)
13+
{
14+
Assert.AreEqual(matrix1.GetLength(0), matrix2.GetLength(0));
15+
Assert.AreEqual(matrix1.GetLength(1), matrix2.GetLength(1));
16+
for (int i = 0; i < matrix1.GetLength(0); i++)
17+
{
18+
for (int j = 0; j < matrix1.GetLength(1); j++)
19+
{
20+
Assert.AreEqual(matrix1[i, j], matrix2[i, j], epsilon, $"At index ({i}, {j})");
21+
}
22+
}
23+
}
24+
25+
public double[,] GenerateRandomMatrix(int m, int n)
26+
{
27+
double[,] result = new double[m, n];
28+
Random random = new Random();
29+
for (int i = 0; i < m; i++)
30+
{
31+
for (int j = 0; j < n; j++)
32+
{
33+
result[i, j] = random.NextDouble() - 0.5;
34+
}
35+
}
36+
return result;
37+
}
38+
39+
[Test]
40+
public void RandomUnitVector()
41+
{
42+
double epsilon = 0.0001;
43+
// unit vector should have length 1
44+
Assert.AreEqual(1, V.Magnitude(ThinSvd.RandomUnitVector(10)), epsilon);
45+
// unit vector with single element should be [-1] or [+1]
46+
Assert.AreEqual(1, Math.Abs(ThinSvd.RandomUnitVector(1)[0]), epsilon);
47+
// two randomly generated unit vectors should not be equal
48+
Assert.AreNotEqual(ThinSvd.RandomUnitVector(10), ThinSvd.RandomUnitVector(10));
49+
}
50+
51+
[Test]
52+
public void Svd_Decompose()
53+
{
54+
CheckSvd(new double[,] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } });
55+
CheckSvd(new double[,] { { 1, 2, 3 }, { 4, 5, 6 } });
56+
CheckSvd(new double[,] { { 1, 0, 0, 0, 2 }, { 0, 3, 0, 0, 0 }, { 0, 0, 0, 0, 0 }, { 0, 2, 0, 0, 0 } });
57+
}
58+
59+
[Test]
60+
public void Svd_Random([Random(3, 10, 5)] int m, [Random(3, 10, 5)] int n)
61+
{
62+
double[,] matrix = GenerateRandomMatrix(m, n);
63+
CheckSvd(matrix);
64+
}
65+
66+
public void CheckSvd(double[,] testMatrix)
67+
{
68+
double epsilon = 1E-5;
69+
double[,] u;
70+
double[,] v;
71+
double[] s;
72+
(u, s, v) = ThinSvd.Decompose(testMatrix, 1E-8, 1000);
73+
74+
for (int i = 1; i < s.Length; i++)
75+
{
76+
// singular values should be arranged from greatest to smallest
77+
Assert.GreaterOrEqual(s[i - 1], s[i]);
78+
}
79+
80+
for (int i = 0; i < u.GetLength(1); i++)
81+
{
82+
double[] extracted = new double[u.GetLength(0)];
83+
// extract a column of u
84+
for (int j = 0; j < extracted.Length; j++)
85+
{
86+
extracted[j] = u[j, i];
87+
}
88+
89+
if (s[i] > epsilon)
90+
{
91+
// if the singular value is non-zero, then the basis vector in u should be a unit vector
92+
Assert.AreEqual(1, V.Magnitude(extracted), epsilon);
93+
}
94+
else
95+
{
96+
// if the singular value is zero, then the basis vector in u should be zeroed out
97+
Assert.AreEqual(0, V.Magnitude(extracted), epsilon);
98+
}
99+
}
100+
101+
for (int i = 0; i < v.GetLength(1); i++)
102+
{
103+
double[] extracted = new double[v.GetLength(0)];
104+
// extract column of v
105+
for (int j = 0; j < extracted.Length; j++)
106+
{
107+
extracted[j] = v[j, i];
108+
}
109+
110+
if (s[i] > epsilon)
111+
{
112+
// if the singular value is non-zero, then the basis vector in v should be a unit vector
113+
Assert.AreEqual(1, V.Magnitude(extracted), epsilon);
114+
}
115+
else
116+
{
117+
// if the singular value is zero, then the basis vector in v should be zeroed out
118+
Assert.AreEqual(0, V.Magnitude(extracted), epsilon);
119+
}
120+
}
121+
122+
// convert singular values to a diagonal matrix
123+
double[,] expanded = new double[s.Length, s.Length];
124+
for (int i = 0; i < s.Length; i++)
125+
{
126+
expanded[i, i] = s[i];
127+
}
128+
129+
130+
// matrix = U * S * V^t, definition of Singular Vector Decomposition
131+
AssertMatrixEqual(testMatrix, u.Multiply(expanded).Multiply(M.Transpose(v)), epsilon);
132+
AssertMatrixEqual(testMatrix, u.Multiply(expanded.Multiply(M.Transpose(v))), epsilon);
133+
}
134+
}
136135
}

0 commit comments

Comments
 (0)