We address the problem of identifying a graph structure from the observation of signals defined on its nodes. Fundamentally, the unknown graph encodes direct relationships between signal elements, which we aim to recover from observable indirect relationships generated by a diffusion process on the graph. The fresh look advocated here permeates benefits from convex optimization and stationarity of graph signals, in order to identify the graph shift operator (a matrix representation of the graph) given only its eigenvectors. These spectral templates can be obtained, e.g., from the sample covariance of independent graph signals diffused on the sought network. The novel idea is to find a graph shift that, while being consistent with the provided spectral information, endows the network with certain desired properties such as sparsity. To that end we develop efficient inference algorithms stemming from provably-tight convex relaxations of natural nonconvex criteria, particularizing the results for two shifts: the adjacency matrix and the normalized Laplacian. Algorithms and theoretical recovery conditions are developed not only when the templates are perfectly known, but also when the eigenvectors are noisy or when only a subset of them are given. Numerical tests showcase the effectiveness of the proposed algorithms in recovering social, brain, and amino-acid networks.