This article presents a new hybrid method (combining search based methods and direct construction methods) to generate all $4 \times 4$ involutory maximum distance separable (MDS) matrices over $\mathbf{F}_{2^m}$. The proposed method reduces the search space complexity at the level of $$\sqrt n $$, where n represents the number of all $4 \times 4$ invertible matrices over $\mathbf{F}_{2^m}$ to be searched for. Hence, this enables us to generate all $4 \times 4$ involutory MDS matrices over $\mathbf{F}_{2^3}$ and $\mathbf{F}_{2^4}$. After applying global optimization technique that supports higher Exclusive-OR (XOR) gates (e.g., XOR3, XOR4) to the generated matrices, to the best of our knowledge, we generate the lightest involutory/non-involutory MDS matrices known over $\mathbf{F}_{2^3}$, $\mathbf{F}_{2^4}$ and $\mathbf{F}_{2^8}$ in terms of XOR count. In this context, we present new $4 \times 4$ involutory MDS matrices over $\mathbf{F}_{2^3}$, $\mathbf{F}_{2^4}$ and $\mathbf{F}_{2^8}$, which can be implemented by 13 XOR operations with depth 5, 25 XOR operations with depth 5 and 42 XOR operations with depth 4, respectively. Finally, we denote a new property of Hadamard matrix, i.e., (involutory and MDS) Hadamard matrix form is, in fact, a representative matrix form that can be used to generate a small subset of all $2^k\times 2^k$ involutory MDS matrices, where k > 1. For k = 1, Hadamard matrix form can be used to generate all involutory MDS matrices.