Generalized linear models (GLMs) are fundamental tools for statistical modeling, with maximum likelihood estimation (MLE) serving as the classical approach for parameter inference. While MLE performs well for canonical GLMs, it can become computationally challenging in more general settings with non-canonical, non-smooth, or nonlinear link functions, where the resulting optimization landscape may be ill-conditioned, non-convex, or non-differentiable. In this paper, we study an alternative estimation framework based on variational inequalities (VIs), which formulates GLM estimation through an operator-based equilibrium condition rather than likelihood minimization. We analyze the VI estimator from a statistical perspective and establish finite-sample error bounds and asymptotic normality under mild regularity conditions, together with convergence guarantees for fixed-point and stochastic approximation algorithms. The framework accommodates a broad class of link functions, including non-canonical and non-monotone cases satisfying a strong Minty-type condition, and extends naturally to generalized additive models via basis expansion. Numerical experiments demonstrate that the VI approach achieves competitive finite-sample accuracy and improved numerical stability relative to MLE, particularly in GLMs and GAMs with non-canonical or non-smooth link functions.