Numerical optimization algorithms often require the (symmetric) matrix of second derivatives, $\nabla^{2} f (x)$, of some problem function $f: R^{n} \rightarrow R^{1}$. If the Hessian matrix is large and sparse then estimation by finite differences can be quite attractive since several schemes allow for estimation in $\ll n$ gradient evaluations. The purpose of this paper is to analyze, from a combinatorial point of view, a class of methods known as substitution methods. We present a concise characterization of such methods in graph-theoretic terms. Using this characterization, we develop a complexity analysis of the general problem and derive a roundoff error bound on the Hessian approximation. Moreover, the graph model immediately reveals procedures to effect the substitution process optimally (ie. using fewest possible substitutions given the differencing directions) in space proportional to the number of nonzeroes in the Hessian matrix.