There's a simple way by which _"each round of DES algorithm is its own inverse"_. Consider round $n$ of DES as involving (almost  only) a function $g_n$ with
$$g_n(L\mathbin\|R)=\bigl(L\oplus f(R,K_n)\bigr)\mathbin\|R$$
where $K_n$ is the 48-bit subkey for round $n$, function $f$ is the _"cipher function"_ (given in the [definition of DES][1]), and $L$ and $R$ are 32-bit bitsrings forming a 64-bit block.

That function $g_n(L\mathbin\|R)$ verifies $g_n(g_n(L\mathbin\|R))=L\mathbin\|R$, as thought in the question; or in other words $g_n$ is an [involution][2]; or in yet other words $g_n\circ g_n$ is the identity function. Proof:
$$\begin{align}
g_n(g_n(L\mathbin\|R))&=g_n\Bigl(\bigl(L\oplus f(R,K_n)\bigr)\mathbin\|R\Bigr)\\
&=\Bigl(\bigl(L\oplus f(R,K_n)\bigr)\oplus f(R,K_n)\Bigr)\mathbin\|R\\
&=\Bigl(L\oplus\bigl(f(R,K_n)\oplus f(R,K_n)\bigr)\Bigr)\mathbin\|R\\
&=\left(L\oplus0^{32}\right)\mathbin\|R\\
&=L\mathbin\|R
\end{align}$$
That proof invokes the definition of $g_n$ (twice), [associativity][3] of $\oplus$, that $f$ is a function, that for all 32-bit $X$ it holds $X\oplus X=0^{32}$ (the bitstring of 32 zero bits), which is the neutral for $\oplus$.
___
DES encryption chains these 33 operations on 64-bit quantities:
$$\mathsf{IP}\,,\,g_1\,,\,S\,,\,g_2\,,\,S\,,\,\ldots\,,\,S\,,\,g_{15}\,,\,S\,,\,g_{16}\,,\,\mathsf{IP}^{-1}$$
where $S$ is the "swap" involution defined by $S(L\mathbin\|R)=R\mathbin\|L$, function $\mathsf{IP}$ is some permutation of bits, and $\mathsf{IP}^{-1}$ is the inverse permutation.

DES decryption chains these 33 operations on 64-bit quantities:
$$\mathsf{IP}\,,\,g_{16}\,,\,S\,,\,g_{15}\,,\,S\,,\,\ldots\,,\,S\,,\,g_2\,,\,S\,,\,g_1\,,\,\mathsf{IP}^{-1}$$

We see that DES encryption then decryption is the identity function: the $(34-j)^\text{th}$ operation of decryption cancels the $j^\text{th}$ operation of encryption:

- For $j=33$, because $\mathsf{IP}$ cancels $\mathsf{IP}^{-1}$.
- For $j=1$, because $\mathsf{IP}^{-1}$ cancels $\mathsf{IP}$.
- For other even $j$, because $g_{j/2}$ is an involution.
- For other (odd) $j$, because $S$ is an involution.

Importantly, encryption and decryption use the very same structure, only the indexes (that, is, the order of the subkeys $K_n$) differ. That allow to use identical hardware or code for both encryption and decryption.
___
The usual definition of a round of a Feistel cipher typically includes the swap $S$:
$$\begin{align}
g'_n(L\mathbin\|R)&=S\bigl(g_n(L\mathbin\|R)\bigr)\\
&=R\mathbin\|\bigl(L\oplus f(R,K_n)\bigr)
\end{align}$$
and this function is _not_ its own inverse per the sense in the question.

With that notation, DES encryption and decryption are the 18 operations
$$\mathsf{IP}\,,\,g'_1\,,\,g'_2\,,\,\ldots\,,\,g'_{15}\,,\,g_{16}\,,\,\mathsf{IP}^{-1}\\
\mathsf{IP}\,,\,g'_{16}\,,\,g'_{15}\,,\,\ldots\,,\,g'_2\,,\,g_1\,,\,\mathsf{IP}^{-1}$$
In that presentation, it is less apparent that decryption undoes encryption. But that still holds, and becomes obvious again if we expand the $g'_n$ into $g_n$ followed by $S$.

When it is said that DES has 16 rounds, it is ignored $\mathsf{IP}$ and $\mathsf{IP}^{-1}$. And, for a strict [Vulcan][4], only the 16<sup>th</sup> round (the one without swap) is its own inverse.
___
The [DES specification][5], and many textbooks on Feistel ciphers, give indexes to $L$ and $R$ before and after round $n$ of encryption, with $L_0\mathbin\|R_0$ the plaintext after $\mathsf{IP}$. The above equation for $g'_n$ can then be written as:
$$\begin{align}
L_n&\gets R_{n-1}\\
R_n&\gets L_{n-1}\oplus f(R_{n-1},K_n)
\end{align}$$

DES (which as $m=16$ rounds) and some textbooks consider the ciphertext (before $\mathsf{IP}^{-1}$) to be $L_m\mathbin\|R_m$, and specialize the last encryption round's equations:
$$\begin{align}
L_m&\gets L_{m-1}\oplus f(R_{m-1},K_m)\\
R_m&\gets R_{m-1}
\end{align}$$

<sup>There are other conventions around: some texts use the same equations for all rounds. Then some add a final swap, or consider the ciphertext to be $R_m\mathbin\|L_m$; others consider $L_m\mathbin\|R_m$ as the ciphertext (that later kind does not obtain the same ciphertext as DES).</sup>
___
As to decryption, DES uses the same round equations for encryption and decryption except for the numbering of subkeys, that is for $1\le n<m$:
$$\begin{align}
L_n&\gets R_{n-1}&&&L_m&\gets L\oplus f(R_{m-1},K_1)\\
R_n&\gets L\oplus f(R_{n-1},K_{16-n})&&&R_m&\gets R_{m-1}
\end{align}$$
with the decryption's $L_0\mathbin\|R_0$ defined as the encryption's $L_n\mathbin\|R_n$, from which it follows that for $1\le n<m$, the decryption's $L_n\mathbin\|R_n$ is the encryption's $R_{m-n}\mathbin\|L_{m-n}$ (notice the inversion), and the decryption's $L_m\mathbin\|R_m$ is the encryption's $L_0\mathbin\|R_0$.

<sup>Other texts use the same naming for equal variables in decryption and encryption, and define different round equations for decryption. For regular round structure that gives:
$$\begin{align}
R_{n-1}&\gets L_n\\
L_{n-1}&\gets R_n\oplus f(L_n,K_n)
\end{align}$$</sup>

  [1]: https://csrc.nist.gov/csrc/media/publications/fips/46/3/archive/1999-10-25/documents/fips46-3.pdf#page=18
  [2]: https://en.wikipedia.org/wiki/Involution_(mathematics)
  [3]: https://en.wikipedia.org/wiki/Associative_property
  [4]: https://en.wikipedia.org/wiki/Vulcan_(Star_Trek)
  [5]: https://csrc.nist.gov/csrc/media/publications/fips/46/3/archive/1999-10-25/documents/fips46-3.pdf#page=14