Given the ground state wavefunction for an interacting lattice model, we define a correlation density matrix(CDM) for two disjoint, separated clusters $A$ and $B$, to be the density matrix of their union, minus the direct product of their respective density matrices. The CDM can be decomposed systematically by a numerical singular value decomposition, to provide a systematic and unbiased way to identify the operator(s) dominating the correlations, even unexpected ones.