A classical model for the covariance of a random matrix is the Kronecker product of two smaller covariance matrices associated with the rows and columns. Maximum likelihood estimation in such structures involves a non-convex optimization problem and is traditionally handled via an alternating maximization Flip-Flop technique. We prove that the problem is actually geodesically convex on the manifold of positive definite matrices. This allows for better analysis, efficient numerical solutions, and simple extensions through the use of additional geodesically convex regularizations. As an example, we propose to shrink the solution towards the identity when the number of samples is insufficient. We demonstrate the advantages of this approach using computer simulations.