Is there a mechanism of dealing with infinite matrices in Mathematica, and in particular, with linear operators, corresponding to those matrices?
For instance, automatic simplification of operator expressions: $e^D-1 \to \text{DifferenceDelta}$, etc?
Simply speaking: No. Infinite matrices are a word only. In practice they "represent" a class of bounded linear maps between two spaces with a basis fixed in both spaces.
In order to make these kind of maps a useful "algebra" over the Complexes, that is a vector space with concatenation of maps @ as multiplication, the two spaces have to be identical, such that the matrices are transfinite limits of square matrices. But these limits represent a set of measure zero only inside the so called algebra B(V), the set of all bounded (B) automorphisms $V \to V$.
Of course, what one has in mind naively, is the very nice set B(H) where H is the standard Hilbert space $l_2$, the vector space of all complex infinite vectors with finite 2-norm. But even here $B(l_2)$ is much to complex as to allow for the set of simple matrix calculation rules.
Really interesting things happen for matrix operators not in B(H). Standard examples are the matrix operators of multiplying by index $I: f_i \to x f(x)$ and the difference operator $D: f_i \to f_i-f_{i-1}$
They are unbounded operators and therefore have a domain of definition given by the condition, that I f and D f are square integrable.
In Fourier analysis of unbounded operators on the standard Hilbert space of square integrable functions in an interval {a,b}, I and D represent $d: f(k) \to k f(k)$ and $x: f(k)\to f'(k)$ as the Fourier images of $x: f(x)\to x f(x)$ and $d: f(x)\to f'(x)$. This fact means, that most formal algebraic compositions of operators do not make any sense.
Thats the boon and bane of naive reckoning of modern mathematical physics. Sometimes it is working (Heisenberg-Schrödinger-von Neumann) and sometimes it produces an infinite staircase of non-converging error corrections (Feynman-Stückelberg-Wightman)
Mathematica has advanced symbolic manipulation capabilities along with many built-in functionality for mathematical operations. I am almost certain there is no built-in functionality specific to general linear operators, but you can use the symbolic capabilities to do what you want. For example, you can use something such as
expr /. E^D -> 1 + DifferenceDelta
or something similar to manipulate linear operators symbolically. You can implement Umbral calculus in a simple way. If you have specific linear operator manipulations you are unsure how to do, then please ask in another question.
