If you have a starting point that satisfies the LMI constraints, it suffices, in theory, just to keep the determinants positive (a change in sign means an eigenvalue became negative). In practice this could fail if a step size was too large and you jump into an region of the matrix space where, say, two eigenvalues became negative.

Is convex relaxation where you take a nap inside an ellipse?

Any symbolic effort in this direction will, at best, only handle toy problems. It's just hugely worse in terms of computational complexity.

One way, probably not the best, would be to check the eigenvalues of the matrix difference.

matrixGreaterEqual[m1_, m2_] /; 
  MatrixQ[m1, Element[#, Reals] &] && 
   MatrixQ[m2, Element[#, Reals] &] && 
   Dimensions[m1] == Dimensions[m2] && Apply[Equal, Dimensions[m1]] :=
  Module[{diffmat = m1 - m2, evals},
  evals = Eigenvalues[diffmat];
  VectorQ[evals, Element[#, Reals] &] && Min[evals] >= 0]

Example:

m1 = {{5, -1}, {-1, 3}};
m2 = {{2, -1}, {-1, 2}};

In[11]:= matrixGreaterEqual[N@m1, m2]

(* Out[11]= True *)

There is probably a better way to go about this using CholeskyDecomposition.

There are matrix power series, a matrix exponential (can be defined in terms of power series), and more. But I should point out that all this is somewhat removed from the intent of this forum.