(Too long for a comment.) I managed to numerically extract the condition numbers up to $n = 30$. I have plotted this on a log scale: [![enter image description here][1]][1] It looks a bit faster than linear on a log scale; maybe a small quadratic term? But it looks no faster than $e^{n^2}$. Code: ```cpp #include <Eigen/Dense> #include <Eigen/Eigenvalues> #include <boost/math/special_functions/factorials.hpp> #include <boost/multiprecision/cpp_bin_float.hpp> #include <iostream> #include <vector> #include <fstream> using Real = boost::multiprecision::cpp_bin_float_100; using Mat = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>; using boost::math::factorial; int main() { std::ofstream of{"data.csv"}; std::cout << std::setprecision(std::numeric_limits<Real>::digits10); of << "n, cond(M)\n"; for (int n = 2; n < 200; ++n) { Mat X = Mat::Random(n,n); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { Real d3 = 2*n - i - j - 1; X(i,j) = 1/(factorial<Real>(n - i - 1)*factorial<Real>(n - j - 1)*d3); } } Eigen::SelfAdjointEigenSolver<Mat> es(n); es.compute(X); std::vector<Real> v(es.eigenvalues().data(), es.eigenvalues().data() + n); if (!std::is_sorted(v.begin(), v.end())) { std::cerr << "Expected invariant is broken.\n"; return 1; } Real cond = v.back()/v.front(); if (cond > 1/std::numeric_limits<Real>::epsilon()) { std::cerr << "Precision must be increased to get more samples; fails at n = " << n << "\n"; break; } of << n << ", " << cond << "\n"; } of.close(); } ``` Generated data: ``` n, cond(M) 2, 19.2815 3, 1181.56 4, 165823 5, 4.18166e+07 6, 1.65669e+10 7, 9.47936e+12 8, 7.39574e+15 9, 7.54511e+18 10, 9.7498e+21 11, 1.55626e+25 12, 3.00702e+28 13, 6.91676e+31 14, 1.86767e+35 15, 5.84992e+38 16, 2.10375e+42 17, 8.60899e+45 18, 3.97753e+49 19, 2.06044e+53 20, 1.18933e+57 21, 7.60721e+60 22, 5.36477e+64 23, 4.15244e+68 24, 3.51305e+72 25, 3.2363e+76 26, 3.23507e+80 27, 3.49783e+84 28, 4.07854e+88 29, 5.11458e+92 30, 6.87997e+96 ``` [1]: https://i.sstatic.net/gosWQ.png