[X,fX,trace_min,flat,error_flat,norm_H,range]=NCtraceOptRand(f,S,d) computes lower bounds for trace of f over all symmetric matrices from D_S: ie bounds L_{\Theta^2}^{s} for s=d/2 and s=d/2+delta, delta=deg(S)/2
Function also tries to extend the optimal solution underlying L_{\Theta^2}^{d/2+delta} to a flat solution using randomization idea

f is an NCpoly representing a polynomial
S is a set of nc polynomails defining D_S
d is a starting degree for the hierarchy (even number)

X: from GNS - a matrix where each of its rows represents a square matrix
fX: f(X) where X is from GNS
trace_min:  L_{\Theta^2}^{s}
flat = 1 if program finds flat extension
flat = 0, if program does not find flat ext.
flat = -1 if the primal problem is infeasible (f is not in the module
error_flat ... norm of the difference between the flat extension returned by randomization idea and the brute force flat extension (which is no more feasible for constraints)
norm_H ... Frobenious norm of flat extension returned by randomization idea
range ... diference between ranks of original matrix and flat extension used 3 methods for rank computation

NCtraceOptRand(f,S,d)