Wrong compute of the matrix rank!


(Lili) #1

For a matrix of 4-by-4, as follows
Matrix A = DenseMatrix.OfArray(new double[,] {
{4,4,1,3},
{1,-2,1,0},
{4,0,2,2},
{7,6,2,5}});
The rank computed with the equation as follow is 3
var ranka = A.Rank();
However, the true rank of the matrix is 2.


#2

Yes, the result is incorrect.
To understand why, need to look at how to calculate the rank.

http://numerics.mathdotnet.com/api/MathNet.Numerics.LinearAlgebra.Double/DenseMatrix.htm#Rank

“effective numerical rank, obtained from SVD”

Going to the github:

public override int Rank
{
   get
   {
      double tolerance = Precision.DoublePrecision*Math.Max(U.RowCount, VT.RowCount);
      return S.Count(t => Math.Abs(t) > tolerance);
   }
}

Let’s look at the SVD:

        Matrix A = DenseMatrix
            .OfArray(new double[,] { 
            { 4, 4, 1, 3 }, 
            { 1,-2, 1, 0 }, 
            { 4, 0, 2, 2 }, 
            { 7, 6, 2, 5 } });
        
        Console.WriteLine(A);
        var svd = A.Svd();
        Console.WriteLine(svd.S);

Print:

DenseMatrix 4x4-Double
4   4  1  3
1  -2  1  0
4   0  2  2
7   6  2  5

DenseVector 4-Double
     13.103
    3.78295
6.26902E-16
          0

https://dotnetfiddle.net/ksmcIf

So, problem in tolerance:

double tolerance = Precision.DoublePrecision*Math.Max(U.RowCount, VT.RowCount);

I tried to find where got this formula and found that:

Take:

s = svd(A);
tol = max(size(A))*eps(max(s));
r = sum(s > tol);

Tried reproduce in Scilab:

So, %eps = Precision.DoublePrecision != eps .

So what is this Eps in Matlab?

http://www.mathworks.com/help/matlab/ref/eps.html

"d = eps(x), where x has data type single or double, returns the positive distance from abs(x) to the next larger floating-point number of the same precision as x. If x has type duration, then eps(x) returns the next larger duration value. The command eps(1.0) is equivalent to eps.

Now can find the answer:

Putting it all together, we obtain:

 static double Eps(double value)
        {
            long bits = BitConverter.DoubleToInt64Bits(value);
            double nextValue = BitConverter.Int64BitsToDouble(bits + 1);
            return nextValue - value;
        }
 static int MyRank(Matrix m)
        {
            var svd = m.Svd();
            double tolerance = Eps(svd.S.Max()) * Math.Max(svd.U.RowCount, svd.VT.RowCount);
            return svd.S.Count(x => Math.Abs(x) > tolerance);
        }
public static void Main(string[] args)
        {
            Matrix A = DenseMatrix
                .OfArray(new double[,] { 
                { 4, 4, 1, 3 }, 
                { 1,-2, 1, 0 }, 
                { 4, 0, 2, 2 }, 
                { 7, 6, 2, 5 } });

            Console.WriteLine(MyRank(A));

        }

https://dotnetfiddle.net/UfE9ME

Print: 2

Sorry for my English.


(Lili) #3

Thank you very much for your kind help and response to my question. I understand and will adopt your good advice.
Ps: you English is very good and easy to understand.

Best Regards.

Lili.

usedsometimes@163.com


(Christoph Rüegg) #4

Thanks for the thorough analysis!

I’ve changed the implementation accordingly to use Precision.EpsilonOf instead of the 1-normalized Precision.DoublePrecision in mainline, (to be) released as part of Math.NET Numerics v3.10.