I think I inadvertently answered your question in your other forum thread. 
To summarize, for y = x^2, the decoders d are solved in such a way such that \sum_{i} a_i d_i approximates the function x^2.
More precisely, the decoders d = \Gamma^{-1}\Upsilon, where \Gamma = \mathbf{A}^T \cdot \mathbf{A} (where \mathbf{A} is the matrix containing the a_i(\mathbf{x}) for all sample points \mathbf{x} and for all neurons i) and \Upsilon = f(\mathbf{x}) \cdot \mathbf{A} (i.e., the dot product of the function f applied to all sample points \mathbf{x} with the matrix \mathbf{A}). For the case of y = x^2, f(\mathbf{x}) = x^2. This is on page 68 of the NEF book.