why not Fitting GLMs with least squares? My question is why we don't use least square to fit Generalized linear model parameters and instead always use maximum likelihood.
 A: GLMs ARE usually fit using iteratively reweighted least squares, see here and references list there, and this post ! This method is based on maximizing the maximum likelihood objective based on Fisher scoring, which is a variant of Newton-Raphson. In a single step one could only approximate the true ML function using least squares though - this would then come down to using a single step of this Fisher scoring algorithm. In practice, the solution obtained in that way can still be quite good though, especially if you use sensible initialisations (e.g. using 1/(y+1) as approximate 1/variance weights in the weighted least squares regression if you are dealing with Poisson noise where the mean=variance). See this post for an example.
A minimal implementation would be:
glm.irls = function(X, y, family=binomial, maxit=25, tol=1e-08, beta.start=rep(0,ncol(X))) {
    if (is.function(family)) family <- family()
    beta = beta.start
    for(j in 1:maxit)
    {
      eta    = X %*% beta
      g      = family$linkinv(eta)
      gprime = family$mu.eta(eta)
      z      = eta + (y - g) / gprime
      W      = as.vector(gprime^2 / family$variance(g))
      betaold   = beta
      beta      = as.matrix(coef(lm.wfit(x=X, y=z, w=W)),ncol=1) # regular weighted LS fit = solve(crossprod(X,W*X), crossprod(X,W*z))
      if(sqrt(crossprod(beta-betaold)) < tol) break
    }
    list(coefficients=beta, iterations=j, z=z, weights=W, X=X, y=y, wlmfit=lm(z~1+X[,-1], weights=W))
}

If you use a distribution with an identity link you can see that in the algorithm above z=y and each iteration just comes down to doing a weighted least squares regression with 1/variance weights. For Poisson e.g. one would then use initial weights = 1/(y+small epsilon) and iterate this based on the predicted yhat, where your weights will then become 1/yhat. With Gaussian errors (ie regular OLS regression) the weights will all just be equal to 1 which means you will get your solution in 1 iteration as there are mo weights to optimize.
Example for logistic regression:
data("Contraception",package="mlmRev")
R_GLM = glm(use ~ age + I(age^2) + urban + livch,
            family=binomial, x=T, data=Contraception)
IRLS_GLM = glm.irls(X=R_GLM$x, y=R_GLM$y, family=binomial)
print(data.frame(R_GLM=coef(R_GLM), IRLS_GLM=coef(IRLS_GLM))) # coefficients match with glm output
                   R_GLM     IRLS_GLM
(Intercept) -0.949952124 -0.949952124
age          0.004583726  0.004583726
I(age^2)    -0.004286455 -0.004286455
urbanY       0.768097459  0.768097459
livch1       0.783112821  0.783112821
livch2       0.854904050  0.854904050
livch3+      0.806025052  0.806025052

A: I can think of a couple reasons: one theoretical, the other practical.
The theoretical reason is because we want to use maximum likelihood estimation (MLE.) This not only gives us a principled way to pose the model fitting procedure as an optimization problem, there are also theoretical results which only work when we use MLE, such as Wilk's theorem which gives rise to the likelihood ratio test and the analysis of deviance, or Wald's test which gives us a way to test the significance of parameters in a generalized linear model. All of these tools go away if we don't use maximum likelihood estimation. 
However, there is also another very practical answer: the squared error loss function is not always convex for GLMs. For example, it is easy to prove this true for logistic regression. You can even see the loss function start to curve downward near 0 and 1:

So even if we wanted to use MSE as an atheoretical loss function this non-convexity could cause problems during fitting!
