mathnet-numerics icon indicating copy to clipboard operation
mathnet-numerics copied to clipboard

Why is the determinant of the zero matrix NaN?

Open neomaru opened this issue 8 years ago • 5 comments

Matrix m1=new DenseMatrix(3,3,new double[]{ 0,0,0, 0,0,0, 0,0,0 }); Console.WriteLine(m1); Console.WriteLine("det={0}",m1.Determinant()); //det=NaN

//OR

Matrix m2=new DenseMatrix(3,3,new double[]{ 0,0,0, 1,1,1, 1,1,1 }); Console.WriteLine(m2); Console.WriteLine("det={0}",m2.Determinant()); //det=NaN

//However, the following result is correct

Matrix m3=new DenseMatrix(3,3,new double[]{ 0.5,0,0, 0 ,1,0, 0 ,0,1 }); Console.WriteLine(m3); Console.WriteLine("det={0}",m3.Determinant());//det=0.5

neomaru avatar Dec 01 '16 06:12 neomaru

It's actually not happening. Take a look at the print screen below image

gandarez avatar Dec 08 '16 15:12 gandarez

Hello gandarez,

Thank you for your verification.

Later, I also checked and found the phenomena seems to happen when MKL(x86 or x64) library is used. This error does not happen in managed library. This problem does not happen in MKL library when the determinant is not zero.

My execution environment is as follows:

Windows10 1607 VisualStudio 2015 Proffessional Update3 Target Project .NET 4.5 (The same problem happens both in virtual environment and Windows7.)

I added following from nuget package manager with console application.

MathNet.Numerics v3.13.1 MathNet.Numerics.MKL.Win-x64 v2.2.0

Would you ask your colleague who uses Windows environment to verify?

The code I used is as follows.

using System;
using System.Collections.Generic;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;

namespace ConsoleApplication1
{
    class Program
    {
        static void Main(string[] args)
        {
            Control.UseManaged();
            Console.WriteLine(Control.LinearAlgebraProvider.ToString());

            Matrix<double> m1 = new DenseMatrix(3, 3, new double[]{
                0,0,0,
                0,0,0,
                0,0,0
            });
            Console.WriteLine("Managed det={0}", m1.Determinant()); // 'Managed det=0' correct

            Control.UseNativeMKL();
            Console.WriteLine(Control.LinearAlgebraProvider.ToString());

            Matrix<double> m2 = new DenseMatrix(3, 3, new double[]{
                0,0,0,
                0,0,0,
                0,0,0
            });
            Console.WriteLine("MKL det={0}", m2.Determinant()); // 'MKL det=NaN' wrong!!

            Matrix<double> m3 = new DenseMatrix(3, 3, new double[]{
                0.5,0,0,
                0,  1,0,
                0,  0,1
            });
            Console.WriteLine("MKL det={0}", m3.Determinant()); // 'MKL det=0.5' When the determinant becomes nonzero, the result is correct.
        }
    }
}

math net_error

neomaru avatar Dec 09 '16 14:12 neomaru

I figured out the problem, the vector becomes NaN when calling an external method d_lu_factor at line 276 on MathNet.Numerics.Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider.Double.cs

I haven't figured out yet where MKL library is generated and how it works since there's no managed code.

I'm not sure if the correct approach is to return 0 when gets NaN value.

public override double Determinant
        {
            get
            {
                var det = 1.0;
                for (var j = 0; j < Factors.RowCount; j++)
                {
                    if (Pivots[j] != j)
                    {
                        det *= -Factors.At(j, j);
                    }
                    else
                    {
                        det *= Factors.At(j, j);
                    }
                }

                return double.IsNaN(det) ? 0 : det;
            }
        }

gandarez avatar Dec 09 '16 18:12 gandarez

Maybe http://www.netlib.org/lapack/lapacke.html#_nan_checking

gandarez avatar Dec 09 '16 20:12 gandarez

I am against replacing a NaN result simply with a 0. This is because in some cases the result really should be NaN: when one of the values in the matrix is NaN (think: unknown) then the result must be NaN (unknown) too.

There are only two proper ways to fix this:

  • file a bug report for MKL and wait till they have fixed it
    • drawback: can take a long time and might not be fixed for backward compatibility
  • check the input for all zeros and if so, return zero, otherwise let MKL calculate it
    • drawback: checking for all zeros has a performance impact if you have lots/large sparse matrices

MovGP0 avatar Dec 25 '16 04:12 MovGP0