Not a new solution but an elaboration from a computational viewpoint.

From a practical point of view one can solve your problem as follows:

Firstly, we do not care about the diagonal values being $1$.
Secondly, we enlarge the problem to real symmetric matrices with
no strictly negative eigenvalues.

Such a matrix $A$ of size $p\times p$ is then given by $A=\sum_i \lambda_i v_i^tv_i$
where $v_1,\dots,v_p$ is an orthogonal basis of eigenvectors with eigenvalues $\lambda_1\geq \lambda_2\geq\dots$
of $A$.

Suppose there is a fast way for computing an eigenvector $v_1$ (of norm $1$)
associated to the largest eigenvalue $\lambda_1$ of $A$. (One can for instance consider the projective
limit of $A^k v$ for $v\in\mathbb R^p$ a generic vector.)

The largest eigenvalue $\lambda_1$ is then given by
$$\lambda_1=\sum_{i,j}A_{i,j}{P(1)}_{i,j}=\sum_{i,j}A_{i,j}{(v_1)}_i{(v_1)}_j$$
where $P(1)=v_1^t v_1$ is the orthogonal rank one projector associated to $v_1$.

Replacing $A$ by $A-\lambda_1P_1$ and iterating one gets $P_2,P_3,\dots$ and $\lambda_2,\lambda_3,\dots$.

Your solution is then given by $\sum_{i=1}^p\lambda_iP(i)=\sum_{i=1}^p\lambda_iv_i^tv_i$.