Condition number for a symmetric positive definite matrix

Question

I am trying to get an estimate for the induced 2-norm condition number $\kappa_{2}(M)$ of this matrix $M$: $$M_{ij} = \frac{1}{(n-i)!(n-j)!(2n-i-j+1)} = \displaystyle\int_{0}^{1}\frac{x^{n-i}}{(n-i)!} \frac{x^{n-j}}{(n-j)!}\:{\rm{d}}x,$$ which is symmetric positive definite of size $n\times n$. The above link gives the elements of the inverse. We also know that its determinant is the reciprocal of this OEIS sequence.

I want to obtain something like $\kappa_{2}(M) = \Omega(f(n))$ for some to-be-determined $f(n)$.

From numerics, it seems that $\kappa_2$ grows very fast with respect to $n$. But I am not sure what might be an analytic way to obtain the asymptotic for $\lambda_{\max}/\lambda_{\min}$.

Answer 1

(Too long for a comment.)

I managed to numerically extract the condition numbers up to $n = 112$. I have plotted this on a log scale:

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}$. The best fit line is dashed. The best fit parabola of the natural logarithm of the condition number is $-33.5838 + 8.09351n + 0.0217842n^2$. I also tried fitting a cubic polynomial as a sanity check, obtaining $-19.6223 + 6.71268n + 0.0516761n^2 - 0.000174806n^3$; that the coefficient of $n^3$ is negative and small might be indicative that it isn't "real", but of course we can't rule that out numerically.

Code:

#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 <cmath>
#include <fstream>

using Real = boost::multiprecision::number<boost::multiprecision::cpp_bin_float<1000>>;
using Mat = Eigen::Matrix<Real, Eigen::Dynamic, Eigen::Dynamic>;
using boost::math::factorial;
using std::log;
using std::sqrt;

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) {
        std::cout << "n = " << n << "\n";
        Mat X = Mat::Zero(n,n);
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                Real d3 = Real(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;
        }
        assert(v.front() > 0);

        Real cond = v.back()/v.front();
        if (cond > 0.01/sqrt(std::numeric_limits<Real>::epsilon())) {
            std::cerr << "Precision must be increased to get more samples; fails at n = " << n << "\n";
            break;
        }
        of << n << ", " << cond << std::endl;
    }
    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.88024e+96
31, 9.90493e+100
32, 1.52261e+105
33, 2.49406e+109
34, 4.34471e+113
35, 8.03433e+117
36, 1.57442e+122
37, 3.26413e+126
38, 7.14851e+130
39, 1.65132e+135
40, 4.01803e+139
41, 1.02846e+144
42, 2.76577e+148
43, 7.80505e+152
44, 2.30875e+157
45, 7.15071e+161
46, 2.31656e+166
47, 7.84207e+170
48, 2.7714e+175
49, 1.02155e+180
50, 3.92398e+184
51, 1.56944e+189
52, 6.5308e+193
53, 2.82524e+198
54, 1.26967e+203
55, 5.92328e+207
56, 2.86664e+212
57, 1.43826e+217
58, 7.4761e+221
59, 4.02365e+226
60, 2.24086e+231
61, 1.29065e+236
62, 7.68353e+240
63, 4.72539e+245
64, 3.00063e+250
65, 1.96638e+255
66, 1.3292e+260
67, 9.26348e+264
68, 6.65305e+269
69, 4.92194e+274
70, 3.74915e+279
71, 2.93919e+284
72, 2.37052e+289
73, 1.96612e+294
74, 1.67632e+299
75, 1.46866e+304
76, 1.32174e+309
77, 1.22145e+314
78, 1.15867e+319
79, 1.12785e+324
80, 1.12618e+329
81, 1.15317e+334
82, 1.21051e+339
83, 1.30226e+344
84, 1.43536e+349
85, 1.62041e+354
86, 1.87313e+359
87, 2.21652e+364
88, 2.68421e+369
89, 3.32574e+374
90, 4.21478e+379
91, 5.46219e+384
92, 7.23698e+389
93, 9.80034e+394
94, 1.35617e+400
95, 1.91725e+405
96, 2.76844e+410
97, 4.08212e+415
98, 6.14524e+420
99, 9.44279e+425
100, 1.48074e+431
101, 2.36913e+436
102, 3.86672e+441
103, 6.43655e+446
104, 1.09254e+452
105, 1.89067e+457
106, 3.33506e+462
107, 5.99552e+467
108, 1.09826e+473
109, 2.04958e+478
110, 3.89608e+483
111, 7.54267e+488
112, 1.4869e+494

(These number break every plotting program I've tried past $n = 75$ as that's where they start overflowing double.)