This paper introduces a causality constrained sequential matrix diagonalisation (SMD) algorithm, which generates a causal paraunitary transformation that approximately diagonalises and spectrally majorises a parahermitian matrix, and can be used to determine a polynomial eigenvalue decomposition. This algorithm builds on a multiple shift technique which speeds up diagonalisation per iteration step based on a particular search space, which is constrained to permit a maximum number of causal time shifts. The results presented in this paper show the performance in comparison to existing algorithms, in particular an unconstrained multiple shift SMD algorithm, from which our proposed method derives.