Generalized methods of moments (GMM) with many indicator variables and ill-conditioned covariance matrix of the error terms

I'm having a problem estimating the optimal weighting matrix for the second stage of (feasible efficient) two-step GMM. Let me explain what I am doing:

For the first stage I use a consistent but inefficient estimate of the optimal weighting matrix (W) (1/N(Z'Z))^-1, where Z is a N x L matrix of included and excluded instruments (L) for N observations and estimate

beta_hat = (X'Z W Z'X)^-1*(X'Z W Z'y)

where X is a matrix of independent variables (K) with dimension N x K and y is the dependent variable.

To get the optimal weighting matrix I calculate the residual from the first stage GMM estimates u_hat= y-X*beta_hat and form the block-diagonal matrix Omega_hat = (u_hat_1^2 ... u_hat_N^2). The optimal weighting matrix is then given by W_hat = (1/N(Z' Omega_hat Z))^-1. ( Where 1/N(Z' Omega_hat Z) is the covariance matrix of the error terms)

Now then I try to calculate this with Matlab. I often get the error

W_opt = inv((Z'*omega_hat*Z)/N)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.921073e-19.

What could the reason be for this behavior. Might it be problematic that my matrix Z contains about 200+ indicator variables which are often zero for the 40000+ observations in my data set? Apart from these indicator variables there are only about 20 "normal" included and excluded instruments in Z.

Now I have read this blog post by Loren Shure on what inv() is for and she notes that for inv(X'X) we should use QR decomposition instead [Purpose of inv]. So why code actually currently looks more like this:

u2diag= spdiags(u_hat.*u_hat,0,length(u),length(u));
sqrtu2diagZ= (u.^.5)*IV;
[~,R] = qr(sqrtu2diagZ,0);
S = inv(R);
W_opt = S*S'*N;

But this does not always help.

So my questions is, should I expect that 1/N(Z' Omega_hat Z) is ill-conditioned given my many indicator variables?