Closed form inverse for the matrix $\exp(-|m-n|)$ In the context of physically coupled oscillators with short-range interactions, the following interaction potential appears:
$V(x_1,x_2) = e^{-a|x_1-x_2|}$
Here, $x_1,x_2 \in \mathbb R$ and $a>0$.
In calculating some of the physical properties of the system, the expression $$\left[\delta(x_1-x_2)-b V(x_1,x_2)\right]^{-1}$$ needs to be calculated, where $b\neq-1$. Here $\delta$ is the Dirac-delta distribution. This effectively requires knowledge of the functional inverse of $V$, or $V^{-1}(x_1,x_2)$
I have stepped back from this problem and attempted to first examine the presumably simpler case of a discretized system.
Now $x_1,x_2 \in \{-L,-L+1,...,L-1, L\}$. Then I have the square matrix $$V_{x_1,x_2} = e^{-a|x_1-x_2|}$$
Or perhaps more simply:
$$V_{m,n} = e^{-a|m-n|}$$
This matrix is real, positive-definite, symmetric, and Toeplitz so it is definitely invertible and seemingly very well-behaved. For large enough values of $a$ the matrix approaches the identity matrix exponentially fast.
The question I have is whether this special matrix has a closed form solution for its inverse. Or, at least, whether this special matrix has ever been studied in any detail. By eye, the matrix looks so symmetric I feel it must have appeared in another context.
Ex. For an 4x4 system
$$V=\begin{pmatrix}
 1&  e^{-a}&  e^{-2a}& e^{-3a}\\ 
 e^{-a}&  1&  e^{-a}& e^{-2a}\\ 
 e^{-2a}&  e^{-a}&  1& e^{-a}\\ 
 e^{-3a}&  e^{-2a}&  e^{-a}& 1
\end{pmatrix}$$
Edit: I attempted to invert the above matrix with Mathematica to get the following matrix.
$$V=\frac{-e^{a}}{1-e^{2a}}\begin{pmatrix}
 e^{a}&  -1&  0& 0\\ 
 -1&  e^{a}+e^{-a}&  -1& 0\\ 
 0&  -1&  e^{a}+e^{-a}& -1\\ 
 0&  0&  -1& e^{a}
\end{pmatrix}$$
 A: Let $\eta$ be the $n \times n$ matrices with superdiagonal entries $1$ and $0$ elsewhere. Let $\bar{\eta}$ be its transpose. 
$$\eta_{ij} = \begin{cases} 1, & j = i + 1\\0, &\text{otherwise}\end{cases}
\quad\text{ and }\quad 
\bar{\eta}_{ij} = \begin{cases} 1, & j = i - 1\\0, &\text{otherwise}\end{cases}
$$
They are nilpotent matrices ($\eta^n = \bar{\eta}^n = 0$ ) and satisfy identities
$$\eta \bar{\eta} \eta = \eta \quad\text{ and }\quad\bar{\eta} \eta \bar{\eta} = \bar{\eta}\tag{*1}$$
It is easy to verify both $\bar{\eta}\eta$ and $\eta\bar{\eta}$ are diagonal matrices with all but one diagonal entries equal to $1$.
$\bar{\eta}\eta$ vanishes on the first diagonal entry while $\eta\bar{\eta}$ vanishes on the last diagonal entry.
$$(\bar{\eta}\eta)_{ij} = \begin{cases} 1, & i = j > 1\\0, & \text{otherwise}
\end{cases}
\quad\text{ and }\quad
(\eta\bar{\eta})_{ij} = \begin{cases} 1, & i = j < n\\0, & \text{otherwise}\end{cases}
$$
Let $\lambda = e^{-a}$. In terms of $\eta$ and $\bar{\eta}$, the matrix $V$ at hand can be rewritten as
$$V = I + \sum_{k=1}^{n-1} \lambda^k \eta^k  + \sum_{k=1}^{n-1}\lambda^k \bar{\eta}^k
= I + \sum_{k=1}^{\infty} \lambda^k \eta^k  + \sum_{k=1}^{\infty}\lambda^k \bar{\eta}^k
= (I - \lambda \eta)^{-1} + (I - \lambda\bar{\eta})^{-1}\lambda\bar{\eta}$$
Multiply $V$ by $(I - \lambda\bar{\eta})$ on the left and $(I-\lambda\eta)$ on the right, we get
$$(I - \lambda\bar{\eta})V(I-\lambda\eta) = (I - \lambda\bar{\eta}) + \lambda\bar{\eta}(I - \lambda\eta) = I - \lambda^2\bar{\eta}\eta\tag{*2}$$
The rightmost matrix is a diagonal matrix with $1$ at first
diagonal entry and $1 - \lambda^2$ at the rest of the diagonal.
When $\lambda^2 \ne 1$, it is invertible with inverse 
$$\left(I - \lambda^2 \bar{\eta}\eta\right)^{-1} = I + \frac{\lambda^2}{1-\lambda^2}\bar{\eta}\eta\tag{*3}$$
Invert both sides of $(*2)$ and apply $(*1)$ and $(*3)$, we obtain
$$V^{-1} = ( 1 - \lambda\eta)\left( I + \frac{\lambda^2}{1-\lambda^2}\bar{\eta}\eta \right)( 1 - \lambda \bar{\eta})
= I -\frac{\lambda}{1-\lambda^2}(\eta + \bar{\eta}) + \frac{\lambda^2}{1-\lambda^2}\left(\bar{\eta}\eta + \eta\bar{\eta}\right)
$$
This means $V^{-1}$ is a tridiagonal matrix with 
$\frac{1}{1-\lambda^2}$ at first and last diagonal entries,
$\frac{1+\lambda^2}{1-\lambda^2}$ at the remaining diagonal entries
and $-\frac{\lambda}{1-\lambda^2}$ on the subdiagonals. i.e $V^{-1}$ has the form
$$
V^{-1} = \frac{1}{1-\lambda^2}\begin{bmatrix}
1 & -\lambda & 0 &  0 & \ldots & 0 & 0 & 0 & 0 \\
-\lambda & 1 + \lambda^2 & -\lambda & 0 & \ldots & 0 & 0 & 0 & 0\\
0 & -\lambda & 1 + \lambda^2 & -\lambda & \ldots & 0 & 0 & 0 & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0 & \ldots & -\lambda & 1+\lambda^2 & -\lambda & 0\\
0 & 0 & 0 & 0 & \ldots & 0 & -\lambda & 1+\lambda^2 & -\lambda\\
0 & 0 & 0 & 0 & \ldots & 0 & 0 & - \lambda & 1 
\end{bmatrix}
$$
A: Your computation for $V^{-1}$ already suggests that $V^{-1}$ behaves like a Laplacian on your discrete lattice. Let us check this by the following heuristic computation:
\begin{align*}
\frac{d}{dx} V(x, y) &= -a\operatorname{sign}(x - y)V(x,y), \\
\frac{d^2}{dx^2} V(x, y) &= -2a\delta(x-y) + a^2V(x,y).
\end{align*}
This suggests that $ V^{-1} = \frac{1}{2a}\left( a^2 - \frac{d^2}{dx^2} \right)$. To make this argument a bit more rigorous, we investigate $V$ on the Fourier side. Notice first that
$$ \int_{-\infty}^{\infty} \frac{2a}{a^2+4\pi^2 t^2} \, e^{2\pi i x t} \, dt = e^{-a|x|}. $$
Thus for $\varphi \in L^1(\Bbb{R})$, Fubini's theorem tells that
\begin{align*}
V\varphi(x)
&= \int_{-\infty}^{\infty} V(x, y)\varphi(y) \, dy \\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{2a}{a^2+4\pi^2 t^2} \, e^{2\pi i (x-y) t}\varphi(y) \, dt \, dy \\
&= \int_{-\infty}^{\infty} \underbrace{\frac{2a}{a^2+4\pi^2 t^2}}_{=:\hat{V}(t)} \, e^{2\pi i x t} \mathcal{F}\varphi(t) \, dt
 = \mathcal{F}^{-1}(\hat{V}(\mathcal{F}\varphi))(x)
\end{align*}
That is, $V$ is simply a multiplicative operator on the Fourier side. So its inverse is given by
$$ V^{-1}\varphi(x) = \int_{-\infty}^{\infty} \frac{a^2+4\pi^2 t^2}{2a} \, e^{2\pi i x t} \mathcal{F}\varphi(t) \, dt, $$
which we easily identify as the operator $\frac{1}{2a}\left(a^2 - \frac{d^2}{dx^2}\right)$ as claimed.

For the discrete case, Mathematica suggests that
$$ [V^{-1}]_{m,n} = \frac{1}{e^{a} - e^{-a}} \times \begin{cases}
e^a,         & i = j \in \{-L, L\} \\
e^a + e^{-a} & i = j \in \{-L+1, \cdots, L-1\} \\
-1           & |i - j| = 1 \\
0            & |i - j| \geq 2
\end{cases} $$
which I have no simple idea to prove at this point. Notice that this guess is again consistent with our previous computation in the following sense: Consider the operator $V(x,y)$ on the rescaled lattice $\{ -L\delta, \cdots, L\delta \}$, then with a nice function $\varphi$, we have
\begin{align*}
&\frac{1}{\delta}\sum_{n} [V^{-1}]_{m \delta, n \delta} \cdot \varphi(n \delta) \\
&\hspace{2em} = \frac{1}{\delta\sinh(a\delta)}\left( \cosh(a\delta)\varphi(mx) - \frac{\varphi((m-1)\delta) + \varphi((m+1)\delta)}{2} \right) \\
&\hspace{4em} \xrightarrow[\quad \substack{m\delta &\to x \\ \delta &\to 0} \quad]{} \frac{1}{2a}(a^2\varphi(x) - \varphi''(x)).
\end{align*}
