LA: Matrix Moore-Penrose Pseudo-Inverse (SVD) #432

src/Numerics/LinearAlgebra/Complex/Matrix.cs

@@ -450,7 +450,45 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             Map(Complex.Log, result, Zeros.Include);
         }
 
+        /// <summary>
+        /// Computes the Moore-Penrose Pseudo-Inverse of this matrix.
+        /// </summary>
+        public override Matrix<Complex> PseudoInverse()
+        {
+            var svd = Svd(true);
+            var w = svd.W;
+            var s = svd.S;
+            double tolerance = Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.DoublePrecision;
 
+            for (int i = 0; i < s.Count; i++)
+            {
+                s[i] = s[i].Magnitude < tolerance ? 0 : 1/s[i];
+            }
+
+            w.SetDiagonal(s);
+            return (svd.U * w * svd.VT).Transpose();
+        }
+
+        /// <summary>
+        /// Computes the trace of this matrix.
+        /// </summary>
+        /// <returns>The trace of this matrix</returns>
+        /// <exception cref="ArgumentException">If the matrix is not square</exception>
+        public override Complex Trace()
+        {
+            if (RowCount != ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
+            }
+
+            var sum = Complex.Zero;
+            for (var i = 0; i < RowCount; i++)
+            {
+                sum += At(i, i);
+            }
+
+            return sum;
+        }
 
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
@@ -638,27 +676,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return Vector<Complex>.Build.Dense(ret);
         }
 
-        /// <summary>
-        /// Computes the trace of this matrix.
-        /// </summary>
-        /// <returns>The trace of this matrix</returns>
-        /// <exception cref="ArgumentException">If the matrix is not square</exception>
-        public override Complex Trace()
-        {
-            if (RowCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSquare);
-            }
-
-            var sum = Complex.Zero;
-            for (var i = 0; i < RowCount; i++)
-            {
-                sum += At(i, i);
-            }
-
-            return sum;
-        }
-
         /// <summary>
         /// Evaluates whether this matrix is hermitian (conjugate symmetric).
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/Matrix.cs

@@ -445,6 +445,25 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             Map(Complex32.Log, result, Zeros.Include);
         }
 
+        /// <summary>
+        /// Computes the Moore-Penrose Pseudo-Inverse of this matrix.
+        /// </summary>
+        public override Matrix<Complex32> PseudoInverse()
+        {
+            var svd = Svd(true);
+            var w = svd.W;
+            var s = svd.S;
+            float tolerance = (float)(Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.SinglePrecision);
+
+            for (int i = 0; i < s.Count; i++)
+            {
+                s[i] = s[i].Magnitude < tolerance ? 0 : 1/s[i];
+            }
+
+            w.SetDiagonal(s);
+            return (svd.U * w * svd.VT).Transpose();
+        }
+
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()

src/Numerics/LinearAlgebra/Double/Matrix.cs

@@ -416,6 +416,46 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             Map(Math.Log, result, Zeros.Include);
         }
 
+        /// <summary>
+        /// Computes the Moore-Penrose Pseudo-Inverse of this matrix.
+        /// </summary>
+        public override Matrix<double> PseudoInverse()
+        {
+            var svd = Svd(true);
+            var w = svd.W;
+            var s = svd.S;
+            double tolerance = Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.DoublePrecision;
+
+            for (int i = 0; i < s.Count; i++)
+            {
+                s[i] = s[i] < tolerance ? 0 : 1/s[i];
+            }
+
+            w.SetDiagonal(s);
+            return (svd.U * w * svd.VT).Transpose();
+        }
+
+        /// <summary>
+        /// Computes the trace of this matrix.
+        /// </summary>
+        /// <returns>The trace of this matrix</returns>
+        /// <exception cref="ArgumentException">If the matrix is not square</exception>
+        public override double Trace()
+        {
+            if (RowCount != ColumnCount)
+            {
+                throw new ArgumentException(Resources.ArgumentMatrixSquare);
+            }
+
+            var sum = 0.0;
+            for (var i = 0; i < RowCount; i++)
+            {
+                sum += At(i, i);
+            }
+
+            return sum;
+        }
+
         /// <summary>Calculates the induced L1 norm of this matrix.</summary>
         /// <returns>The maximum absolute column sum of the matrix.</returns>
         public override double L1Norm()
@@ -602,27 +642,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return Vector<double>.Build.Dense(ret);
         }
 
-        /// <summary>
-        /// Computes the trace of this matrix.
-        /// </summary>
-        /// <returns>The trace of this matrix</returns>
-        /// <exception cref="ArgumentException">If the matrix is not square</exception>
-        public override double Trace()
-        {
-            if (RowCount != ColumnCount)
-            {
-                throw new ArgumentException(Resources.ArgumentMatrixSquare);
-            }
-
-            var sum = 0.0;
-            for (var i = 0; i < RowCount; i++)
-            {
-                sum += At(i, i);
-            }
-
-            return sum;
-        }
-
         /// <summary>
         /// Evaluates whether this matrix is hermitian (conjugate symmetric).
         /// </summary>

src/Numerics/LinearAlgebra/Matrix.Arithmetic.cs

@@ -1629,6 +1629,9 @@ namespace MathNet.Numerics.LinearAlgebra
             return LU().Inverse();
         }
 
+        /// <summary>Computes the Moore-Penrose Pseudo-Inverse of this matrix.</summary>
+        public abstract Matrix<T> PseudoInverse();
+
         /// <summary>
         /// Computes the Kronecker product of this matrix with the given matrix. The new matrix is M-by-N
         /// with M = this.Rows * lower.Rows and N = this.Columns * lower.Columns.

src/Numerics/LinearAlgebra/Single/Matrix.cs

@@ -416,6 +416,25 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             Map(x => (float)Math.Log(x), result, Zeros.Include);
         }
 
+        /// <summary>
+        /// Computes the Moore-Penrose Pseudo-Inverse of this matrix.
+        /// </summary>
+        public override Matrix<float> PseudoInverse()
+        {
+            var svd = Svd(true);
+            var w = svd.W;
+            var s = svd.S;
+            float tolerance = (float)(Math.Max(RowCount, ColumnCount) * svd.L2Norm * Precision.SinglePrecision);
+
+            for (int i = 0; i < s.Count; i++)
+            {
+                s[i] = s[i] < tolerance ? 0 : 1/s[i];
+            }
+
+            w.SetDiagonal(s);
+            return (svd.U * w * svd.VT).Transpose();
+        }
+
         /// <summary>
         /// Computes the trace of this matrix.
         /// </summary>