In this paper, we introduce two novel split least-squares mixed element procedures for pseudo-parabolic equations. By selecting the least-squares functional properly, each procedure can be split into two independent symmetric positive definite sub-procedures. One of sub-procedures is for the primitive unknown variable u, which is the same as the standard Galerkin finite element procedure and the other is for the introduced flux variable σ. Optimal order error estimates are developed. A numerical example is given to show the efficiency of the introduced schemes.