We consider the problem of finding a low-rank approximate solution to a system of linear equations in symmetric, positive semidefinite matrices. Specifically, let $A_1,\ldots,A_m \in \R^{n\times n}$ be symmetric, positive semidefinite matrices, and let $b_1,\ldots,b_m \ge 0$. We show that if there exists a symmetric, positive semidefinite matrix $X$ to the system $A_i \bullet X = b_i$ (where $i = 1,\ldots,m$), then for any fixed $d = 1,\ldots,O(\log m)$, there exists an $X_0 \succeq 0$ of rank at most $d$ such that $\beta \cdot b_i \le A_i \bullet X_0 \le \alpha \cdot b_i$ for all $i = 1,\ldots,m$, where $\alpha = 1+O(\log m/d)$; $\beta = \Omega(m^{-2/d})$ for $d = O(\log m / \log\log m)$, and $\beta = \Omega((\log m)^{-3\log m / (d\log\log m)})$. Moreover, such an $X_0$ can be found in randomized polynomial time. This complements a result of Barvinok and provides a unified treatment of and generalizes several results in the literature.