Does this generalization of matrix exponential formula make sense? In a linear algebra class we defined the exponential of some invertible linear matrix $M$ through a power series
$$
\exp M =\sum_{k=0}^\infty \frac{M^k}{k!},
$$
where if one had to write the explicit matrix indices they would have
$$
[\exp M]_{ij} =\text{I}_{ij}+M_{ij}+\frac{1}{2!} \sum_m M_{im}M_{mj}+\frac{1}{3!} \sum_{mn} M_{im}M_{mn}M_{nj}+\cdots,
$$
here the discrete sums run over the dimension of the matrix $M$.
I was wondering if there is a continuum generalization of such sum/expression, where the discrete sum is replaced by integrals over some domain $D$. Say that $M$ becomes some function of two variables $m(x, y)$, then what is the sum (and when does it converge?)
$$
\delta(x-y)+m(x,y)+\frac{1}{2!}\int_D dz\, m(x,z)m(z,y)
+\frac{1}{3!}\int_{D^2} dz\, dw\, m(x,z)m(z,w)m(w, y)+\cdots$$
where in place of the identity matrix I wrote Dirac's delta function.
Is this completely crazy or not? Does the sum equal some sort of $\exp m$ ?
Thank you!
 A: Matrices are discrete objects, and matrix exponential is defined through Taylor series for the exponential function.
Matrices represent linear operators on finite-dimensional vector spaces. I guess the continuous analogue would be to consider linear operators on uncountably infinite-dimensional vector spaces, like spaces of functions. These naturally come with some extra structure - often you can take limits of a series of functions, for instance - and they're studied in functional analysis.
In your question you just replaced a matrix with a function in two variables, and this makes no sense to me (but waiting for other users to correct me!) in virtue of what I wrote above, thinking about your generalisation from a matrix to an operator.
If you don't think in terms of matrices, then you can replace $e^{M}$ with $e^{f(x, y)}$ (or you can replace the whole thing with another two variables funtion $h(x, y)$). In this case you will make use of Taylor series in $2$ variables, when expanding your function.
In your case you have finite sums because you're computing matrix products, which added together will give you the exponential matrix.
If you do the same with a $2$ variables function, you won't have any "index sum" as you do with matrices. All the single terms will add up to the infinite to give your initial function (or it's value at a point, provided the series converges).
If you integrate term by term, you just get the integral of your initial function (under certain conditions).
If it's an exponential then you'll get
$$e^{f(x, y)} \approx e^{f(x_0, y_0)} + e^{f(x_0, y_0)} f'(x_0, y_0) + \ldots$$
For other function you'll obey to the general Taylor series to get the required series.
Passing from discrete con continuum is always a tricky game, maybe not always possible, especially if you just step from a matrix to a $2$ variables function.
For example, why would you change the identity operator with the Heaviside generalised function (aka Dirac Delta distribution)? I know that it's the equivalent in the continuum of the Kronecker $\delta$, but at that point when using Dirac distribution, your sum of integrals may lack of sense, for when $x = y$ it's all good (ish), but when $x\neq y$ you get a singularity.
A: Here's an attempt to formally relate your integral formula to the traditional matrix exponential.
Suppose that $m:[0,1]\times [0,1] \to \Bbb C$. We are interested in the function $f:[0,1]\times[0,1] \to\Bbb C$ given by
$$
f(x,y) = \delta(x-y)+m(x,y)+\frac{1}{2!}\int_D dz_1\, m(x,z_1)m(z,y)
+\cdots
$$
We can get an approximation for this quantity using the value of $f$ over $x,y \in \{0,1/n,\dots,1-1/n\}$ by using Riemann approximations for the integrals. In particular, we have
$$
\frac 1{k!} \int_{D^{k-1}} m(x,z_1)m(z_{k-1},y)\prod_{j=1}^{k-2}m(z_j,z_{j+1}) \,dz_1 \cdots dz_k
\\ \approx \frac 1{k!}n^{1-k} \sum_{0 \leq p_1,\dots,p_{k-1} \leq n-1} m(x,p_1/n) m(p_{k-1}/n,y) \prod_{j=1}^{k-2}m(p_j/n,p_{j+1}/n).
$$
If $M$ has entries $M_{ij} = m(i/n,j/n)$ for $1 \leq i,j \leq n-1$ and we have $x = i/n, y = j/n$ for some $i,j$, then the above approximation is the $i,j$ of the matrix $\frac {1}{k! \,n^{k-1}} M^{k}$. Thus, in order to approximate $f(x,y)$, we could compute the sum
$$
I + M + \frac 1{2!\,n} M^2 + \frac{1}{3!\,n^2}M^3 + \cdots \\
= I + n\left[(M/n) + \frac 1{2!}(M/n)^2 + \frac 1{3!} (M/n)^3 + \cdots\right]\\
= I + n[\exp(M/n) - I].
$$
So, perhaps your function can be interpreted as a "limit" of this matrix as $n \to \infty$.
A: The way you should see this is as follows:

*

*Matrices are linear operators that map (finite-dimensional) vectors onto (finite-dimensional) vectors:
$$
  y = Mx.
$$
As a consequence, you use two indices, one for $x$ and one for $y$, and you can write the formula above as
$$
  y_i = \sum_j M_{ij} x_j.
$$
Computing the matrix exponential $e^M$ then works as you describe above, namely: It is defined via a Taylor series.


*There is no conceptual problem to extending this approach to operators that map from countably-infinite to countably-infinite spaces -- think vectors with infinitely many entries. Examples of such spaces are $\ell_2$, or functions that can be written as Fourier series -- for the latter, you can represent the function by its (infinitely many) Fourier coefficients, and operators that map functions in such a space to another function in that space are simply infinite-dimensional matrices for which the matrix exponential is defined in the same way. The only difference is that the sums run over infinitely many values of $i, j, k$.


*Your question is then what happens when you go to operators that map from and to spaces that do not have a countable basis. Say, for example, that you are considering an operator $\cal M$ that maps one function $f\in L^\infty$ to another function $g\in L^\infty$:
$$
  g = {\cal M} f,
$$
in the sense that if you evaluate $g$ at a point $x$, you have
$$
  g(x) = ({\cal M} f)(x),
$$
which is not (necessarily) the same as
$$
  g(x) = {\cal M} (f(x)),
$$
unless $\cal M$ is local. An example of a non-local operator of this kind would be a convolution integral
$$
({\cal M} f)(x) = \int_D m(x,y) f(y) \, dy,
$$
but it could also be local differentiation:
$$
({\cal M} f)(x) = \frac{df}{dx}(x).
$$
What, then, is $e^{\cal M}$? By Taylor expansion, it is
$$
  e^{\cal M} = I + {\cal M} + \frac 1{2!} {\cal M}^2  + \frac 1{3!} {\cal M}^3 + \cdots.
$$
If the operator can be written as a convolution integral with some $m$, this yields
$$
  (e^{\cal M} f)(x) 
  = 
  If + 
  {\cal M} f
+ \frac 1{2!} {\cal M}^2 f
+ \frac 1{3!} {\cal M}^3 f
  + \cdots
\\
  = 
  \int_D \delta(x-y) f(y) \, dy + 
  \int_D m(x,y) f(y) \, dy
+ \frac 1{2!} \int_{D}\int_D m(x,y) m(y,z) f(z) \, dz \, dy
\\
+ \frac 1{3!} \int_{D}\int_D\int_D m(x,y) m(y,z) m(z,a) f(a) \, da \, dz \, dy
  + \cdots.
$$
As a consequence, the operator $e^{\cal M}$ can be represented as a convolution operator itself, namely
$$
  e^{\cal M}f = \int_D n(x,y) f(y) \, dy
$$
where the "matrix elements" of $e^{\cal M}$ are then
$$
  n(x,y) = 
  \delta(x-y) + 
  m(x,y)
+ \frac 1{2!} \int_D m(x,y) m(y,z) \, dy
\\
+ \frac 1{3!} \int_{D}\int_D m(x,y) m(y,z) m(z,a) \, dz \, dy
  + \cdots.
$$
This is the proper generalization of the matrix exponential. Whether the sum in the last formula converges is something that will depend on the specific properties of $\cal M$, but one can at least write it down.
