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$, and $f$ is the "cipher function" (given in the definition of DES), 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; 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 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 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 $$\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}$$ When it is said that DES has 16 rounds, it is ignored $\mathsf{IP}$ and $\mathsf{IP}^{-1}$ and for a strict Vulcan, only the 16th round (the one without swap) is its own inverse.
The DES specification, 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}$, so that the above equations for $g'_n$ are equivalently worded: $$\begin{align} L_n&\gets R_{n-1}\\ R_n&\gets L\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}$$
But there are other conventions around: some texts use the same equations for all rounds, but in order to get the same result as DES either add a final swap, or consider the ciphertext to be $R_n\mathbin\|L_n$. Yet others use a regular round structure and consider $L_n\mathbin\|R_n$ as the ciphertext.
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$ equals the encryption's $R_{m-n}\mathbin L_{m-n}$ and the decryption's $L_m\mathbin R_m$ equal the encryption's $L_0\mathbin R_0$.
Other texts prefer to have the same naming of equal variables for decryption and encryption, and thus define different round equations for decryption. For a text using a 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}$$