Eigenvectors of the matrix defined by $(i,j) \mapsto (i+j)!$

Question

I have an $n \times n$ integer Hankel matrix $A$ that is defined via

$$A_{i,j} = \left(i+j\right)!$$

where $0!=1$ as is standard and both indices $i$ and $j$ start at $0$. Is there an analytic way to determine the eigenvectors for this Hankel matrix?

The context is that this is relevant to a function approximation problem I'm working on. I wish to use the analytic result so that I can determine an appropriate value of some exponentials in the next step to counterbalance the very large value that will appear, and in doing so make this useful.

---

Edit: The underlying problem was solved, so this specific problem does not need to be solved. Here's the code updated with the polynomials Jair provided:

import numpy as np
def basis_generator(t,order):
    """
    Generate a basis valid over [0,inf), evaluated at points t
Parameters
----------
t : array_like
    Set of positive real points the basis will be generated at. 
order : int
    The maximal order of polynomial to be used in the basis.

Returns
-------
y : array_like
    Orthonormal basis over the points t.
"""
y = np.zeros(shape=(len(t),order+1))
y[:,0] = np.exp(-t/2)
y[:,1] = (1-t)*np.exp(-t/2)
dy = -t*np.exp(-t/2)
for k in range(1,order):
    dy -= dy/(k+1)
    dy -= (t*(y[:,k])/(k+1))
    y[:,k+1] = y[:,k]+dy
return y


#Example 1:
dt = 0.00001
t = np.arange(dt,50,dt)
y = basis_generator(t,8)
print('Condition number: ', np.linalg.cond(y))
print('Inner products: ', dt*(y.T @ y))
#Example 2:
dt = 0.0001
t = np.arange(dt,50,dt)
y = basis_generator(t,100)
f = np.sin(t/5)
b = y.T @ f * dt
print('Relative squared error: ', (np.linalg.norm(f-(y@b))2)/(np.linalg.norm(f)2))

Answer 1

This is not a complete answer, since it does not compute the eigenvectors of the matrix $A$, but it answers a question from the comments.

The Laguerre polynomials $L_n$ defined by

$$L_n(x) = \sum_{k=0}^n \frac{(-1)^k}{k!} {n \choose k} x^k$$

are orthonormal with respect to the measure $e^{-x} d\mu(x)$ on $[0, \infty)$, meaning that $$\int_0^\infty e^{-x} L_i(x) L_j(x) \, dx = \delta_{ij}$$ where $\delta_{ij} = 0$ if $i \neq j$ and $1$ if $i = j$.

Since $\int_0^\infty x^{i+j} e^{-x} dx = (i+j)!$, it follows that

$$L^TAL = I$$ where $L$ is the $(n+1) \times (n+1)$ "Laguerre matrix" of order $n$ with $L_{m,k} = \frac{(-1)^k}{k!} {m \choose k}$ for $0 \leq k \leq m \leq n$ and $A_{ij} = (i+j)!$ for $0 \leq i,j \leq n$. This gives an explicit Cholesky decomposition

$$A = L^{-T}L^{-1}.$$

However, it's not clear to me if it's possible to get an explicit form for the eigenvalues or eigenvectors of $A$.

As an aside, I would note that it would not be hard to find the formula for the orthonormal polynomials with a bit of sleuthing. If you perform Gram-Schmidt on the sequence $1, x, x^2, \ldots$ with respect to the inner product $\langle p,q \rangle = \int_0^\infty p(x) q(x) e^{-x} \, dx$ then you get a solution, and by looking up the coefficients in OEIS you would learn about the Laguerre polynomials.