We optimize the matrix representation of the nucleon-pair approximation (NPA) of the nuclear shell model. The NPA is a widely adopted truncation approach of the nuclear shell model and proves to be effective in describing low-lying states of medium-heavy and heavy nuclei. Due to simplified (yet flexible) commutators and absolute elimination of angular momentum coupling, the matrix representation provides a formalism for the M -scheme NPA more efficient than others as far as we know. It also enables the practicable organization and storage design for intermediate results, including generated collective pairs, matrix products, and matrix traces, so that further optimization is achieved by reducing repetitive matrix operations, which are the most time-consuming procedures in the matrix-represented M -scheme NPA. We also describe optimizations specified for the M -scheme NPA, realized by invoking the Wigner-Eckart theorem, timereversal symmetry, and conjugate operation of spherical tensors. Our optimization makes the combination of matrix representation and NPA more profitable. Such an implementation denoted by optimized matrix representation of NPA (OMR-NPA) is publicly released with open source. Its performance is analyzed and compared against unoptimized NPA codes.