Linear Algebra: Matrix Rank should use effective epsilon #334

Previously it was using the normalized epsilon at 1.0.
See http://discuss.mathdotnet.com/t/wrong-compute-of-the-matrix-rank/120
for discussion.

src/Numerics/LinearAlgebra/Complex/DenseVector.cs

@@ -519,24 +519,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
             return index;
         }
 
-        /// <summary>
-        /// Returns the value of the absolute minimum element.
-        /// </summary>
-        /// <returns>The value of the absolute minimum element.</returns>
-        public override Complex AbsoluteMinimum()
-        {
-            return _values[AbsoluteMinimumIndex()].Magnitude;
-        }
-
-        /// <summary>
-        /// Returns the value of the absolute maximum element.
-        /// </summary>
-        /// <returns>The value of the absolute maximum element.</returns>
-        public override Complex AbsoluteMaximum()
-        {
-            return _values[AbsoluteMaximumIndex()].Magnitude;
-        }
-
         /// <summary>
         /// Returns the index of the absolute maximum element.
         /// </summary>

src/Numerics/LinearAlgebra/Complex/Factorization/Svd.cs

@@ -71,7 +71,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         {
             get
             {
-                double tolerance = Precision.DoublePrecision*Math.Max(U.RowCount, VT.RowCount);
+                double tolerance = Precision.EpsilonOf(S.AbsoluteMaximum().Magnitude)*Math.Max(U.RowCount, VT.RowCount);
                 return S.Count(t => t.Magnitude > tolerance);
             }
         }

src/Numerics/LinearAlgebra/Complex/Vector.cs

@@ -323,7 +323,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex
         /// Returns the value of the absolute minimum element.
         /// </summary>
         /// <returns>The value of the absolute minimum element.</returns>
-        public override Complex AbsoluteMinimum()
+        public sealed override Complex AbsoluteMinimum()
         {
             return At(AbsoluteMinimumIndex()).Magnitude;
         }

src/Numerics/LinearAlgebra/Complex32/DenseVector.cs

@@ -514,24 +514,6 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
             return index;
         }
 
-        /// <summary>
-        /// Returns the value of the absolute minimum element.
-        /// </summary>
-        /// <returns>The value of the absolute minimum element.</returns>
-        public override Complex32 AbsoluteMinimum()
-        {
-            return _values[AbsoluteMinimumIndex()].Magnitude;
-        }
-
-        /// <summary>
-        /// Returns the value of the absolute maximum element.
-        /// </summary>
-        /// <returns>The value of the absolute maximum element.</returns>
-        public override Complex32 AbsoluteMaximum()
-        {
-            return _values[AbsoluteMaximumIndex()].Magnitude;
-        }
-
         /// <summary>
         /// Returns the index of the absolute maximum element.
         /// </summary>

src/Numerics/LinearAlgebra/Complex32/Factorization/Svd.cs

@@ -66,7 +66,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         {
             get
             {
-                double tolerance = Precision.SinglePrecision*Math.Max(U.RowCount, VT.RowCount);
+                double tolerance = Precision.EpsilonOf(S.AbsoluteMaximum().Magnitude)*Math.Max(U.RowCount, VT.RowCount);
                 return S.Count(t => t.Magnitude > tolerance);
             }
         }

src/Numerics/LinearAlgebra/Complex32/Vector.cs

@@ -318,7 +318,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32
         /// Returns the value of the absolute minimum element.
         /// </summary>
         /// <returns>The value of the absolute minimum element.</returns>
-        public override Complex32 AbsoluteMinimum()
+        public sealed override Complex32 AbsoluteMinimum()
         {
             return At(AbsoluteMinimumIndex()).Magnitude;
         }

src/Numerics/LinearAlgebra/Double/DenseVector.cs

@@ -556,24 +556,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double
             return index;
         }
 
-        /// <summary>
-        /// Returns the value of the absolute minimum element.
-        /// </summary>
-        /// <returns>The value of the absolute minimum element.</returns>
-        public override double AbsoluteMinimum()
-        {
-            return Math.Abs(_values[AbsoluteMinimumIndex()]);
-        }
-
-        /// <summary>
-        /// Returns the value of the absolute maximum element.
-        /// </summary>
-        /// <returns>The value of the absolute maximum element.</returns>
-        public override double AbsoluteMaximum()
-        {
-            return Math.Abs(_values[AbsoluteMaximumIndex()]);
-        }
-
         /// <summary>
         /// Returns the index of the absolute maximum element.
         /// </summary>

src/Numerics/LinearAlgebra/Double/Factorization/Svd.cs

@@ -64,7 +64,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         {
             get
             {
-                double tolerance = Precision.DoublePrecision*Math.Max(U.RowCount, VT.RowCount);
+                double tolerance = Precision.EpsilonOf(S.Maximum())*Math.Max(U.RowCount, VT.RowCount);
                 return S.Count(t => Math.Abs(t) > tolerance);
             }
         }

src/Numerics/LinearAlgebra/Single/DenseVector.cs

@@ -546,24 +546,6 @@ namespace MathNet.Numerics.LinearAlgebra.Single
             return index;
         }
 
-        /// <summary>
-        /// Returns the value of the absolute minimum element.
-        /// </summary>
-        /// <returns>The value of the absolute minimum element.</returns>
-        public override float AbsoluteMinimum()
-        {
-            return Math.Abs(_values[AbsoluteMinimumIndex()]);
-        }
-
-        /// <summary>
-        /// Returns the value of the absolute maximum element.
-        /// </summary>
-        /// <returns>The value of the absolute maximum element.</returns>
-        public override float AbsoluteMaximum()
-        {
-            return Math.Abs(_values[AbsoluteMaximumIndex()]);
-        }
-
         /// <summary>
         /// Returns the index of the absolute maximum element.
         /// </summary>

src/Numerics/LinearAlgebra/Single/Factorization/Svd.cs

@@ -64,7 +64,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         {
             get
             {
-                double tolerance = Precision.SinglePrecision*Math.Max(U.RowCount, VT.RowCount);
+                double tolerance = Precision.EpsilonOf(S.Maximum())*Math.Max(U.RowCount, VT.RowCount);
                 return S.Count(t => Math.Abs(t) > tolerance);
             }
         }

src/UnitTests/LinearAlgebraTests/Double/Factorization/SvdTests.cs

@@ -180,6 +180,20 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Double.Factorization
             Assert.AreEqual(factorSvd.Rank, order - 1);
         }
 
+        [Test]
+        public void RankAcceptance()
+        {
+            // http://discuss.mathdotnet.com/t/wrong-compute-of-the-matrix-rank/120
+            Matrix<double> m = DenseMatrix.OfArray(new double[,] {
+                { 4, 4, 1, 3 },
+                { 1,-2, 1, 0 },
+                { 4, 0, 2, 2 },
+                { 7, 6, 2, 5 } });
+
+            Assert.That(m.Svd(true).Rank, Is.EqualTo(2));
+            Assert.That(m.Svd(false).Rank, Is.EqualTo(2));
+        }
+
         /// <summary>
         /// Solve for matrix if vectors are not computed throws <c>InvalidOperationException</c>.
         /// </summary>

src/UnitTests/LinearAlgebraTests/Single/Factorization/SvdTests.cs

@@ -180,6 +180,20 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraTests.Single.Factorization
             Assert.AreEqual(factorSvd.Rank, order - 1);
         }
 
+        [Test]
+        public void RankAcceptance()
+        {
+            // http://discuss.mathdotnet.com/t/wrong-compute-of-the-matrix-rank/120
+            Matrix<float> m = DenseMatrix.OfArray(new float[,] {
+                { 4, 4, 1, 3 },
+                { 1,-2, 1, 0 },
+                { 4, 0, 2, 2 },
+                { 7, 6, 2, 5 } });
+
+            Assert.That(m.Svd(true).Rank, Is.EqualTo(2));
+            Assert.That(m.Svd(false).Rank, Is.EqualTo(2));
+        }
+
         /// <summary>
         /// Solve for matrix if vectors are not computed throws <c>InvalidOperationException</c>.
         /// </summary>