BDCSVD NaN on some 60x60 matrix with EIGEN_USE_BLAS / LLVM / AVX512

Summary

Environment

  • Operating System : Linux
  • Architecture : x86_64
  • Eigen Version : 3.4.0, master
  • Compiler Version : LLVM clang++ 14.0.6 and others LLVM versions
  • Compile Flags : -mavx512f -DEIGEN_USE_BLAS -lblas
  • Vector Extension : AVX512

Minimal Example

eigen_bug.cpp Exerpt:

#include <iostream>
#include "Eigen/Dense"
#include <sstream>

const char* matAsString = R"(
[...] // code with matrix entries is appended
)";

int main()
{
  Eigen::MatrixXd M(60,60);
  std::istringstream iss(matAsString);
  for(int i = 0; i < 60; i++)
    for(int j = 0; j < 60; j++)
      iss >> M(i,j);

  // to check that the matrix was read correctly.
  //std::cout << "M:\n" << M << std::endl;

  auto svd = Eigen::BDCSVD<Eigen::MatrixXd>(M);//, Eigen::ComputeThinV | Eigen::ComputeThinU);
  std::cout << "singular values:\n" << svd.singularValues().transpose() << std::endl;

  return 0;
}

Steps to reproduce

  1. > clang++ -g -Og -mavx512f -DEIGEN_USE_BLAS eigen_bug.cpp -o eigen_bug -lblas
  2. > ./eigen_bug Output:
singular values:
        nan 6.26699e-06 3.04367e-06 6.13745e-08 4.23141e-08 3.62818e-08 1.88088e-08 1.45254e-08 1.08346e-08 7.39115e-09 3.56497e-09 2.48686e-10 1.82483e-10 1.08825e-10 5.89797e-11 4.41218e-11 1.53663e-11           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0

What is the current bug behavior?

The first singular value is NaN, singular values seem to be shifted by one

What is the expected correct behavior?

Correct singular values, no NaNs.

Relevant logs

Correct behavior without AVX512 or with GCC or without EIGEN_USE_BLAS for this specific matrix.

  1. > clang++ eigen_bug.cpp -o eigen_bug
  2. > ./eigen_bug
singular values:
6.26699e-06 3.04367e-06 6.13745e-08 4.23141e-08 3.62818e-08 1.88088e-08 1.45254e-08 1.08346e-08 7.39115e-09 3.56497e-09 2.48686e-10 1.82483e-10 1.08825e-10 5.89797e-11 4.41218e-11 1.53663e-11 1.11324e-20           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0

(I didn't check the resulting singular values, though)

Anything else that might help

I suppose the error is related to almost zero floating point values. I ran this with EIGEN_BDCSVD_DEBUG_VERBOSE and EIGEN_BDCSVD_SANITY_CHECKS enabled:




======================================================================================================================



deflate:9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.62023e-149 1.77998e-150  |  9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.2, set z(1) to zero because 6.09567e-195 < 2.13952e-164  (diag(1)=9.63555e-149)
deflation 4.2, set z(2) to zero because 0 < 2.13952e-164  (diag(2)=9.63555e-149)
deflation 4.2, set z(3) to zero because 0 < 2.13952e-164  (diag(3)=9.63555e-149)
deflation 4.2, set z(5) to zero because 5.28314e-195 < 2.13952e-164  (diag(5)=9.63555e-149)
deflation 4.2, set z(6) to zero because 0 < 2.13952e-164  (diag(6)=9.63555e-149)
deflation 4.2, set z(7) to zero because 4.52706e-180 < 2.13952e-164  (diag(7)=9.63555e-149)
deflation 4.2, set z(8) to zero because 0 < 2.13952e-164  (diag(8)=9.63555e-149)
deflation 4.2, set z(9) to zero because 1.83797e-212 < 2.13952e-164  (diag(9)=9.63555e-149)
deflation 4.2, set z(10) to zero because 0 < 2.13952e-164  (diag(10)=9.63555e-149)
deflation 4.2, set z(11) to zero because 2.88148e-195 < 2.13952e-164  (diag(11)=9.63555e-149)
deflation 4.2, set z(12) to zero because 6.67871e-195 < 2.13952e-164  (diag(12)=9.63555e-149)
deflation 4.2, set z(13) to zero because 1.1625e-179 < 2.13952e-164  (diag(13)=9.62023e-149)
deflation 4.2, set z(14) to zero because 1.46584e-179 < 2.13952e-164  (diag(14)=1.77998e-150)
deflation 4.2, set z(15) to zero because 2.7528e-212 < 2.13952e-164  (diag(15)=9.63555e-149)
deflation 4.2, set z(16) to zero because 0 < 2.13952e-164  (diag(16)=9.63555e-149)
deflation 4.2, set z(17) to zero because 0 < 2.13952e-164  (diag(17)=9.63555e-149)
deflation 4.2, set z(18) to zero because 0 < 2.13952e-164  (diag(18)=9.63555e-149)
deflation 4.2, set z(20) to zero because 0 < 2.13952e-164  (diag(20)=9.63555e-149)
deflation 4.2, set z(21) to zero because 1.66874e-195 < 2.13952e-164  (diag(21)=9.63555e-149)
deflation 4.2, set z(22) to zero because 6.07246e-255 < 2.13952e-164  (diag(22)=9.63555e-149)
deflation 4.2, set z(23) to zero because 0 < 2.13952e-164  (diag(23)=9.63555e-149)
deflation 4.2, set z(24) to zero because 1.31377e-195 < 2.13952e-164  (diag(24)=9.63555e-149)
deflation 4.2, set z(25) to zero because 1.93415e-211 < 2.13952e-164  (diag(25)=9.63555e-149)
deflation 4.2, set z(26) to zero because 1.32706e-211 < 2.13952e-164  (diag(26)=9.63555e-149)
deflation 4.2, set z(27) to zero because 8.99759e-180 < 2.13952e-164  (diag(27)=9.63555e-149)
deflation 4.2, set z(28) to zero because 7.22282e-195 < 2.13952e-164  (diag(28)=9.63555e-149)
to be sorted: 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.62023e-149 1.77998e-150 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149

            :  9.63555e-149             0             0             0 -8.57412e-164             0             0             0             0             0             0             0             0             0             0             0             0             0             0  3.47306e-164             0             0             0             0             0             0             0             0             0

sorted:   [9.6355506e-149, 1.7799783e-150,  9.620229e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149]
      :  9.63555e-149             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0  3.47306e-164             0             0 -8.57412e-164             0             0             0             0             0             0

deflation 4.4 with i = 22 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 21,22 -> 0 -8.57412e-164 0 ; 0 0 -8.57412e-164 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 21 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 20,21 -> 0 0 0 ; 3.47306e-164 0 0 -8.57412e-164
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 20 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 19,20 -> 3.47306e-164 0 0 ; 0 3.47306e-164 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 19 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 18,19 -> 0 3.47306e-164 0 ; 0 0 3.47306e-164 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 18 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 17,18 -> 0 0 0 ; 0 0 0 3.47306e-164
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 16 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 15,16 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 15 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 14,15 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 14 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 13,14 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 13 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 12,13 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 12 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 11,12 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 11 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 10,11 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 10 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 9,10 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 9 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 8,9 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 7 because 9.63555e-149 - 9.63555e-149 == 0 < 2.13952e-164
deflation 4.4: 6,7 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
deflation 4.4 with i = 5 because 9.63555e-149 - 9.63555e-149 == 1.73653e-164 < 2.13952e-164
deflation 4.4: 4,5 -> 0 0 0 ; 0 0 0 0
9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149


j1 =   [9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149,  9.620229e-149, 1.7799783e-150]
j2 =   [9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149, 9.6355506e-149,  9.620229e-149, 1.7799783e-150]

err:      0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# 1

computeSVDofM using:
  z:  9.63555e-149             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0  3.47306e-164             0             0 -8.57412e-164             0             0             0             0             0             0
  d:            0 1.77998e-150 9.62023e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
right-left = 9.63555e-149
     = -1e+12 -99 -24 -10.1111 -5.25 -3.16493 -3 -2.84468 -1.77778 -1.04082 -0.5625 -0.234568 -2e-06
useBisection for k = 0, actual_n = 23
found 9.63555e-149 == 9.63555e-149 + -9.55092e-164 from 0 .. 1.77998e-150
right-left = 0
     = inf inf inf inf inf inf inf inf inf inf inf inf inf
useBisection for k = 19, actual_n = 23
f(0) =inf ; 9.63555e-149 9.63555e-149 9.63555e-149
found 9.63555e-149 == 9.63555e-149 + 1.11254e-308 from 9.63555e-149 .. 9.63555e-149
right-left = 9.63555e-149
     = 2e-06 0.173554 0.305556 0.408284 0.489796 0.54957 0.555556 0.561423 0.609375 0.653979 0.691358 0.722992 0.75
useBisection for k = 22, actual_n = 23
found 9.63555e-149 == 9.63555e-149 + 2.49991e-164 from 9.63555e-149 .. 9.63555e-149
  j:        1.77998e-150 9.62023e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149

  sing-val: 9.63555e-149 1.77998e-150 9.62023e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149
  mu:       -9.55092e-164             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0  1.11254e-308             0             0  2.49991e-164             0             0             0             0             0             0
  shift:    9.63555e-149 1.77998e-150 9.62023e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149 9.63555e-149


    mus:    -9.55092e-164             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0             0  1.11254e-308             0             0  2.49991e-164             0             0             0             0             0             0

    check1 (expect0) : 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

    check2 (>0)      :           1           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0           0 1.80221e-16           0           0           0           0           0           0

zhat(0) =  sqrt( 9.28438e-297)  ;  9.63555e-149 * 9.63555e-149 - 0
     1.98243e-15 == (1.92711e-148 * -9.55092e-164) / (9.63555e-149 * -9.63555e-149)
eigen_bug: ./Eigen/src/SVD/BDCSVD.h:1042: void Eigen::BDCSVD<Eigen::Matrix<double, -1, -1, 0>>::perturbCol0(const Eigen::BDCSVD::ArrayRef &, const Eigen::BDCSVD::ArrayRef &, const Eigen::BDCSVD::IndicesRef &, const Eigen::BDCSVD::VectorType &, const Eigen::BDCSVD::ArrayRef &, const Eigen::BDCSVD::ArrayRef &, Eigen::BDCSVD::ArrayRef) [MatrixType = Eigen::Matrix<double, -1, -1, 0>]: Assertion `prod>=0' failed.
Aborted (core dumped)
Edited by Melven Roehrig-Zoellner