I have a simple matrix equation $$c A + I = (\det A) S$$ which seems linear (or perhaps quadratic) and involves the determinant of the matrix being solved for. Here, $c$ is a constant, and $S$ is a matrix. How could I solve for $A$, a symmetric matrix in $\mathbb{R}^{k\times k}$? An extremely efficient way of computing $A$ would also be fine.
If it helps solve the problem, then know that $S$ corresponds to a (full-rank) covariance matrix $\frac{2}{n}X X^{\intercal}$, and $c$ corresponds to the norm of a mean-vector $\left\Vert \frac{1}{n}\overline{x}\right\Vert_2 = \frac{1}{n^2}\overline{x}^{\intercal}\overline{x}$.
If there is no simple solution, a solution for the $k=2$ case is all I really need.
I'm not sure if this will help or not.
So, let us do case $k=2$.
Let $A$ be a solution and let $Q$ such that $QAQ^{-1}$ is $A$'s Jordan form. Multiplying the equation we get something of form
$$\begin{pmatrix} c\lambda_1+1 & 0\\ 0 & c\lambda_2+1 \end{pmatrix} = \lambda_1 \lambda_2 QSQ^{-1} = \lambda_1 \lambda_2 \begin{pmatrix} s_1 & 0\\ 0 & s_2 \end{pmatrix}$$
which can easily be solved for eigenvalues of $A$. Notice that this works since symmetric matrices are diagonalizable by orthogonal matrix. Also both $A$ and $S$ get into Jordan form under same base change.
Thus, to find solution, diagonalize $S$, calculate eigenvalues of $A$ and return it back: $$A = Q^{-1}\begin{pmatrix} \lambda_1 & 0\\ 0 & \lambda_2 \end{pmatrix}Q.$$
Edit:
To address the issue of finding eigenvalues in higher dimensions.
Since we have $c\lambda_i+1 = s_i\prod\lambda_j$, we can see that $\frac{c\lambda_i+1}{s_i}$ is constant, let us denote it by $d$. So, we have $c\lambda_i+1 =s_id\implies \lambda_i =\frac{ds_i-1}{c}$. Finally, to find $d$, you need to solve $\prod\frac{ds_i-1}{c} = d$, which can be done for $k\leq 4$ by Abel-Ruffini theorem, otherwise you would have to do it numerically.