The loop phase margins of multivariable control systems are defined as the allowable individual loop phase perturbations within which stability of the closed-loop system is guaranteed. This paper presents a frequency domain approach to accurately computing these phase margins for multivariable systems. With the help of unitary mapping between two complex vector space, the MIMO phase margin problem is converted using the Nyquist stability analysis to the problem of some simple constrained optimization, which is then solved numerically with the Lagrange multiplier and Newton-Raphson iteration algorithm. The proposed approach can provide exact margins and thus improves the linear matrix inequalities (LMI) results reported before, which could be conservative.