6
$\begingroup$

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}$.

$\endgroup$
1
  • 2
    $\begingroup$ Numerics seem to suggest that $\lambda_{\text{max}}$ converges to approximately $1.30384272207$ as $n \rightarrow \infty$, so we "just" need to figure out what's happening with $\lambda_{\text{min}}$. $\endgroup$ – Nathaniel Johnston Oct 12 '20 at 1:00
1
$\begingroup$

(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:

enter image description here

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.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.