Molcas developer here, currently implementing symmetric and canonical orthonormalization for the RASSCF module in addition to Gram-Schmidt.
As a rule of thumb you should treat those Orthonormalization schemes (ON-schemes) as different tools for different jobs.
If there are $n$ linear independent atomic orbitals (AOs), of which $n$ molecular orbitals (MOs) are to be constructed, then symmetric ON is probably the best, because the squared deviation of the MOs to the AOs is minimized. This means that your MOs are as close to the AOs as possible, while still being orthonormal.
If there is linear dependency, symmetric orthogonalization can detect, but not cure it.
Canonical and Gram-Schmidt orthonormalization have the advantage, that they can be used to construct only $m$ MOs from $n$ AOs ($n > m$) and to cure linear dependence.
If you compare Gram-Schmidt with Canonical ON, Gram-Schmidt has the disadvantage, that it is not as numerically stable. Besides the AOs that are orthonormalized last, can look completely different after the procedure. If you manually constructed your basis with chemical knowledge, this is not what you want. Modified Gram-Schmidt cures the numerical instability, but still has the disadvantage, that MOs might look completely different to the initial AOs.
PS: Just an addition to the comments: Grahm-Schmidt can be written as linear transformation. Although one would rarely do this.