I asked this question on Mathematica site and the answer was basically "it is impossible". So, some math theory needs to be thrown in besides coding.
Let us represent infinite matrices as formal power series of the form
$$M=\sum _{d=-\infty }^{\infty } m_d(Y)X^d. $$
Here the function $m_d(y)$ represents the element of the matrix $M$ at the diagonal $d$ (positive diagonals are above the main diagonal) with number $y$. Each infinite matrix is uniquely defined by this two-argument function.
In other words, the infinite matrices are basically power series of the following form: $M=\sum_{d=-\infty}^\infty M_d X^d$, where $M_d$ are diaginal matrices and $X=\left( \begin{array}{cccc} 0 & 1 & 0 & . \\ 0 & 0 & 1 & . \\ 0 & 0 & 0 & . \\ . & . & . & . \\ \end{array} \right)$. Since any diafonal matrix can be represented as a function of $Y$: $M_d=m_d(Y)$, where $Y=\left( \begin{array}{cccc} 1 & 0 & 0 & . \\ 0 & 2 & 0 & . \\ 0 & 0 & 3 & . \\ . & . & . & . \\ \end{array} \right)$, we represent them as $M=\sum _{d=-\infty }^{\infty } m_d(Y)X^d$.
This way, for instance, the derivative operator in matrix form would be represented as $D=YX$, integral $\int_0^x f(t)dt$ as $Y^{-1}X^{-1}$, finite difference $\Delta=\exp(D)-1=\sum _{d=1}^{\infty } \frac{(Y X)^d}{d!}=\sum _{d=0}^{\infty } \binom{Y+d-1}{d}X^d-1=\sum _{d=1}^{\infty } \binom{Y+d-1}{d}X^d$.
The symbols $X$ and $Y$ also would represent linear operators: $X:f(x)\to\frac{f(x)}x$ and $Y:f(x)\to xf'(x)$.
That said, I wonder, whether we can define some rules that would show how a function applied to an infinite matrix transforms the kernel function $m_d(Y)$, so that we could manipulate the matrices in the form of such power series in Mathematica?
Am I correct that operations on infinite matrices can be reduced to the operations on formal power series this way?
