Skip to content

Commit 2c63986

Browse files
1fisedi1fisedi
authored
Add Moore-Penrose pseudo-inverse matrix computation algorithm based on SVD (#197)
1 parent 5a6f50e commit 2c63986

6 files changed

Lines changed: 223 additions & 19 deletions

File tree

‎Algorithms.Tests/Numeric/PseudoInverse/PseudoInverseTests.cs‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
using Algorithms.Numeric;
2+
using NUnit.Framework;
3+
using Utilities.Extensions;
4+
5+
namespace Algorithms.Tests.Numeric
6+
{
7+
public static class PseudoInverseTests
8+
{
9+
[Test]
10+
public static void SquaredMatrixInverseWorks()
11+
{
12+
// Arrange
13+
var inMat = new double[,] { { 2, 4, 6 }, { 2, 0, 2 }, { 6, 8, 14 } };
14+
var inMatCopy = new double[,] { { 2, 4, 6 }, { 2, 0, 2 }, { 6, 8, 14 } };
15+
16+
// Act
17+
// using AA+A = A
18+
var result = PseudoInverse.PInv(inMat);
19+
var aainva = inMatCopy.Multiply(result).Multiply(inMatCopy);
20+
21+
var rounded = aainva.RoundToNextInt();
22+
var isequal = rounded.IsEqual(inMatCopy);
23+
// Assert
24+
Assert.IsTrue(isequal);
25+
}
26+
27+
[Test]
28+
public static void NonSquaredMatrixPseudoInverseMatrixWorks()
29+
{
30+
// Arrange
31+
var inMat = new double[,] { { 1, 2, 3, 4 }, { 0, 1, 4, 7 }, { 5, 6, 0, 1 } };
32+
var inMatCopy = new double[,] { { 1, 2, 3, 4 }, { 0, 1, 4, 7 }, { 5, 6, 0, 1 } };
33+
34+
// Act
35+
// using (A+)+ = A
36+
var result = PseudoInverse.PInv(inMat);
37+
var result2 = PseudoInverse.PInv(result);
38+
39+
var rounded = result2.RoundToNextInt();
40+
41+
var isequal = rounded.IsEqual(inMatCopy);
42+
// Assert
43+
Assert.IsTrue(isequal);
44+
}
45+
}
46+
}

‎Algorithms/Numeric/Pseudoinverse/PseudoInverse.cs‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
using System;
2+
using Algorithms.Numeric.Decomposition;
3+
using Utilities.Extensions;
4+
5+
namespace Algorithms.Numeric
6+
{
7+
/// <summary>
8+
/// The Moore–Penrose pseudo-inverse A+ of a matrix A,
9+
/// is a general way to find the solution to the following system of linear equations:
10+
/// ~b = A ~y. ~b e R^m; ~y e R^n; A e Rm×n.
11+
/// There are varios methods for construction the pseudo-inverse.
12+
/// This one is based on Singular Value Decomposition (SVD).
13+
/// </summary>
14+
public static class PseudoInverse
15+
{
16+
/// <summary>
17+
/// Return the pseudoinverse of a matrix based on the Moore-Penrose Algorithm.
18+
/// using Singular Value Decomposition (SVD).
19+
/// </summary>
20+
/// <param name="inMat">Input matrix to find its inverse to.</param>
21+
/// <returns>The inverse matrix approximation of the input matrix.</returns>
22+
public static double[,] PInv(double[,] inMat)
23+
{
24+
// To compute the SVD of the matrix to find Sigma.
25+
var (u, s, v) = ThinSvd.Decompose(inMat);
26+
27+
// To take the reciprocal of each non-zero element on the diagonal.
28+
var len = s.Length;
29+
30+
var sigma = new double[len];
31+
for (var i = 0; i < len; i++)
32+
{
33+
sigma[i] = (Math.Abs(s[i]) < 0.0001) ? 0 : 1 / s[i];
34+
}
35+
36+
// To construct a diagonal matrix based on the vector result.
37+
var diag = sigma.ToDiagonalMatrix();
38+
39+
// To construct the pseudo-inverse using the computed information above.
40+
var matinv = u.Multiply(diag).Multiply(v.Transpose());
41+
42+
// To Transpose the result matrix.
43+
return matinv.Transpose();
44+
}
45+
}
46+
}

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ This repository contains algorithms and data structures implemented in C# for ed
3838
* [Series](./Algorithms/Numeric/Series/)
3939
* [Maclaurin](./Algorithms/Numeric/Series/Maclaurin.cs)
4040
* [Gauss-Jordan Elimination](./Algorithms/Numeric/GaussJordanElimination.cs)
41+
* [Pseudo-Inverse](./Algorithms/Numeric/Pseudoinverse/PseudoInverse.cs)
4142
* [Searches](./Algorithms/Search/)
4243
* [A-Star](./Algorithms/Search/AStar/)
4344
* [Binary](./Algorithms/Search/BinarySearcher.cs)

‎Utilities.Tests/Extensions/MatrixExtensionsTests.cs‎

Lines changed: 60 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -11,12 +11,12 @@ public class MatrixExtensionsTests
1111
public void Multiply_ShouldThrowInvalidOperationException_WhenOperandsAreNotCompatible()
1212
{
1313
// Arrange
14-
var source = new double[,] {{1, 1, 1}, {1, 1, 1}, {1, 1, 1}};
15-
var operand = new double[,] {{1}, {1}};
16-
14+
var source = new double[,] { { 1, 1, 1 }, { 1, 1, 1 }, { 1, 1, 1 } };
15+
var operand = new double[,] { { 1 }, { 1 } };
16+
1717
// Act
1818
Action action = () => source.Multiply(operand);
19-
19+
2020
// Assert
2121
action.Should().Throw<InvalidOperationException>()
2222
.WithMessage("The width of a first operand should match the height of a second.");
@@ -31,7 +31,7 @@ public void Multiply_ShouldCalculateDotProductMultiplicationResult(
3131
public void Copy_ShouldReturnImmutableCopyOfMatrix()
3232
{
3333
// Arrange
34-
var sutMatrix = new double[,] {{1, 1, 1}, {1, 1, 1}, {1, 1, 1}};
34+
var sutMatrix = new double[,] { { 1, 1, 1 }, { 1, 1, 1 }, { 1, 1, 1 } };
3535

3636
// Act
3737
var actualMatrix = sutMatrix.Copy();
@@ -40,42 +40,84 @@ public void Copy_ShouldReturnImmutableCopyOfMatrix()
4040
actualMatrix.Should().NotBeSameAs(sutMatrix);
4141
actualMatrix.Should().BeEquivalentTo(sutMatrix);
4242
}
43-
43+
4444
[Test, TestCaseSource(nameof(MatrixTransposeTestCases))]
4545
public void Transpose_ShouldReturnTransposedMatrix(
4646
double[,] source, double[,] target) =>
4747
source.Transpose().Should().BeEquivalentTo(target);
48-
48+
4949
[Test]
5050
public void MultiplyVector_ShouldCalculateDotProductMultiplicationResult()
5151
{
5252
// Arrange
53-
var source = new double[,] {{2, 2, -1}, {0, -2, -1}, {0, 0, 5}};
54-
var operand = new double[] {2, 2, 3};
55-
var result = new double[] {5, -7, 15};
53+
var source = new double[,] { { 2, 2, -1 }, { 0, -2, -1 }, { 0, 0, 5 } };
54+
var operand = new double[] { 2, 2, 3 };
55+
var result = new double[] { 5, -7, 15 };
5656

5757
// Act
5858
var actualMatrix = source.MultiplyVector(operand);
5959

6060
// Assert
6161
actualMatrix.Should().BeEquivalentTo(result);
6262
}
63-
63+
6464
[Test]
6565
public void Subtract_ShouldThrowArgumentException_WhenOperandsAreNotCompatible()
6666
{
6767
// Arrange
68-
var source = new double[,] {{1, 1, 1}, {1, 1, 1}, {1, 1, 1}};
69-
var operand = new double[,] {{1}, {1}};
70-
68+
var source = new double[,] { { 1, 1, 1 }, { 1, 1, 1 }, { 1, 1, 1 } };
69+
var operand = new double[,] { { 1 }, { 1 } };
70+
7171
// Act
7272
Action action = () => source.Subtract(operand);
73-
73+
7474
// Assert
7575
action.Should().Throw<ArgumentException>()
7676
.WithMessage("Dimensions of matrices must be the same");
7777
}
78-
78+
79+
[Test]
80+
public static void EqualMatricesShouldReturnTrue()
81+
{
82+
// Arrange
83+
var A = new double[,] { { 1, 2, 3 }, { 1, 2, 3 }, { 1, 2, 3 } };
84+
var B = new double[,] { { 1, 2, 3 }, { 1, 2, 3 }, { 1, 2, 3 } };
85+
86+
// Act
87+
var result = A.IsEqual(B);
88+
89+
// Assert
90+
Assert.True(result);
91+
}
92+
93+
[Test]
94+
public static void NonEqualMatricesShouldReturnFalse()
95+
{
96+
// Arrange
97+
var A = new double[,] { { 1, 2, 3 }, { 1, 2, 3 }, { 1, 2, 3 } };
98+
var B = new double[,] { { 1, 2, 3 }, { 1, 2, 6 }, { 1, 2, 3 } };
99+
100+
// Act
101+
var result = A.IsEqual(B);
102+
103+
// Assert
104+
Assert.False(result);
105+
}
106+
107+
[Test]
108+
public static void DifferentSizeMatricesShouldReturnFalse()
109+
{
110+
// Arrange
111+
var A = new double[,] { { 1, 2, 3 }, { 1, 2, 3 }, { 1, 2, 3 } };
112+
var B = new double[,] { { 1, 2, 3 }, { 1, 2, 3 } };
113+
114+
// Act
115+
var result = A.IsEqual(B);
116+
117+
// Assert
118+
Assert.False(result);
119+
}
120+
79121
[Test, TestCaseSource(nameof(MatrixSubtractTestCases))]
80122
public void Subtract_ShouldCalculateSubtractionResult(
81123
double[,] source, double[,] operand, double[,] result) =>
@@ -96,7 +138,7 @@ public void Subtract_ShouldCalculateSubtractionResult(
96138
new double[,] {{11, -22, 29}, {9, -27, 32}, {13, -17, 26}}
97139
}
98140
};
99-
141+
100142
private static readonly object[] MatrixTransposeTestCases =
101143
{
102144
new object[]
@@ -110,7 +152,7 @@ public void Subtract_ShouldCalculateSubtractionResult(
110152
new double[,] {{5, 6}, {8, 9}}
111153
}
112154
};
113-
155+
114156
private static readonly object[] MatrixSubtractTestCases =
115157
{
116158
new object[]

‎Utilities/Extensions/MatrixExtensions.cs‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -128,5 +128,56 @@ public static double[] MultiplyVector(this double[,] matrix, double[] vector)
128128
return result;
129129
}
130130

131+
/// <summary>
132+
/// Perfoms an element by element comparison on both matrices.
133+
/// </summary>
134+
/// <param name="source">Source left matrix.</param>
135+
/// <param name="operand">Openrand right matrix.</param>
136+
/// <returns>true: if all elements are the same; false otherwise.</returns>
137+
public static bool IsEqual(this double[,] source, double[,] operand)
138+
{
139+
if (source.Length != operand.Length ||
140+
source.GetLength(0) != operand.GetLength(0) ||
141+
source.GetLength(1) != operand.GetLength(1))
142+
{
143+
return false;
144+
}
145+
146+
for (var i = 0; i < source.GetLength(0); i++)
147+
{
148+
for (var j = 0; j < source.GetLength(0); j++)
149+
{
150+
if (Math.Abs(source[i, j] - operand[i, j]) >= 0.0001)
151+
{
152+
return false;
153+
}
154+
}
155+
}
156+
157+
return true;
158+
}
159+
160+
/// <summary>
161+
/// Performs a round operation on every element of the input matrix up to the neareast integer.
162+
/// </summary>
163+
/// <param name="source">Input matrix.</param>
164+
/// <returns>Matrix with rounded elements.</returns>
165+
public static double[,] RoundToNextInt(this double[,] source)
166+
{
167+
var rows = source.GetLength(0);
168+
var cols = source.GetLength(1);
169+
170+
var result = new double[rows, cols];
171+
172+
for (var i = 0; i < rows; i++)
173+
{
174+
for (var j = 0; j < cols; j++)
175+
{
176+
result[i, j] = Math.Round(source[i, j]);
177+
}
178+
}
179+
180+
return result;
181+
}
131182
}
132183
}

‎Utilities/Extensions/VectorExtensions.cs‎

Lines changed: 19 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ public static double[] Scale(this double[] vector, double factor)
9191

9292
return result;
9393
}
94-
94+
9595
/// <summary>
9696
/// Transpose 1d row vector to column vector.
9797
/// </summary>
@@ -131,5 +131,23 @@ public static double[] ToRowVector(this double[,] source)
131131

132132
return rowVector;
133133
}
134+
135+
/// <summary>
136+
/// Generates a diagonal matrix from an specified vector.
137+
/// </summary>
138+
/// <param name="vector">The input vector.</param>
139+
/// <returns>A Diagonal matrix</returns>
140+
public static double[,] ToDiagonalMatrix(this double[] vector)
141+
{
142+
var len = vector.Length;
143+
var result = new double[len, len];
144+
145+
for (var i = 0; i < len; i++)
146+
{
147+
result[i, i] = vector[i];
148+
}
149+
150+
return result;
151+
}
134152
}
135153
}

0 commit comments

Comments
 (0)