We describe efficient evaluation of the (global) Hessian matrix of the sum-squared-error measure for CANFIS neuro-fuzzy modular network learning. Our network consists of multiple (local-expert) multilayer perceptrons (MLPs) mediated by fuzzy membership functions, leading to an iteratively reweighted nonlinear least squares problem. In the posed situation, we show how our second-order stagewise backpropagation procedure, recently developed for learning with a single MLP, efficiently exploits the sparsity (of the Hessian matrix) that arises in a multiple-response problem. In spite of its complex modular architecture, our procedure works excellently. Its computational convenience is immense since such an efficient evaluation is crucial in implementing Newton-type second-order algorithms that may exploit negative curvature when the Hessian matrix is indefinite as well as in the Hessian analysis for any type of modular neural-network learning.