In this article, we employ a regression formulation to estimate the high dimensional covariance matrix for a given network structure. Using prior information contained in the network relationships, we model the covariance as a polynomial function of the symmetric adjacency matrix. Accordingly, the problem of estimating a high dimensional covariance matrix is converted to one of estimating low dimensional coefficients of the polynomial regression function, which we can accomplish using ordinary least squares or maximum likelihood. The resulting covariance matrix estimator based on the maximum likelihood approach is guaranteed to be positive definite even in finite samples. Under mild conditions, we obtain the theoretical properties of the resulting estimators. A Bayesian information criterion is also developed to select the order of the polynomial function. Simulation studies and empirical examples illustrate the usefulness of the proposed methods.