Faster matrix completion

[sflammia]

Currently the matrix completion routine uses all the constraints from a maximal Galois orbit and Convex.jl to solve low-rank matrix recovery via solving a semidefinite program. This could probably be made much faster, although it would require care to ensure robustness. Solving via an SDP will not be fast enough if we want to try this in, say, $d = 1000$.

Some possible ideas to speed this up: We should only need about 4d - O(1) constraints for injectivity (this was proven for POVMs by Heinosaari, Mazzarella & Wolf). So instead of using all the constraints from a maximal orbit, we could sample a random sparse set of constraints, say 6d of them to be safe. Since the maximal orbit is generally quite large, sampling might be faster. We could potentially speed this up and remove dependencies on Convex.jl and SCS.jl by using fast non-convex methods like Singular-value thresholding (SVT), Burer-Montiero iterations, FPCS, etc. that come from the theory of low-rank matrix sensing.