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