[update 4] We get an exact result
Here I give a worked example using Pari/GP. I adapted the letters for the matrices, the previously used notation was not very elegant....
We want that
$ \small S = I + A \cdot A^\tau + A^2 \cdot (A^\tau)^2 + A^3 \cdot (A^\tau)^3 + ... $ from where also $ \small S = I + A \cdot S \cdot A^\tau $
The matrix A can be very general; by use of diagonalization of A and the closed form for the geometric series applied to the eigenvalues of A we get an exact result.
[update 5] If A is symmetric, then the process simplifies again: The diagonalization $\small A=W\cdot D \cdot W^{-1} $ gives a matrix W which can be normed to be a rotation-matrix, thus $ \small M=W^{-1}=W^\tau$ , then $ \small Mm = M \cdot M^\tau = I $ , and thus Gg is diagonal with the entries $ \small g_{k,k}= 1/( 1-\lambda_k^2)$ where $\small \lambda_k$ are the eigenvalues (their squares may not be 1 now, but the relation of pairs is now irrelevant). The result is then simply $ \small S = W \cdot Gg \cdot W^\tau $
[/update 5]
\\ Pari/GP d = 4 A = matrix(d,d,r,c,c^r) \\ create some diagonizable matrix W = mateigen(A) M = W^-1 D = diag(M*A*W); mydisplay([W, Mat(D)]) \\ user function for compact output Mm = M*M~ \\ Now the closed-form for the geometric series of the \\ eigenvalues come into play. \\ Note that the eigenvalues( in D) can be larger than 1, \\ but no pair is allowed to give a product of 1 Gg = matrix(d,d, r,c, Mm [r,c] / (1 - D[r] * D[c] ) ) S = W*Gg*W~ \\ exact result \\ then crosscheck by implicite definition chk = matid(d) + A * S * A~ err = round(S - chk,&e) \\ difference is zero
The example-matrices:
$ \qquad \small \begin{array} {lll} A &=& \begin{array} {rrrrr} 1 & 2 & 3 & 4 \\ 1 & 4 & 9 & 16 \\ 1 & 8 & 27 & 64 \\ 1 & 16 & 81 & 256 \end{array} \\ \\ \\ W ,D &=& \begin{array} {rrrr|r} -10.952794 & 8.1645134 & -0.97056263 & 0.017714558 & 0.11272446 \\ 10.508537 & 4.2521778 & -1.8354857 & 0.066926874 & 1.0292480 \\ -5.0996430 & -4.0885199 & -2.6756811 & 0.25725908 & 8.9314971 \\ 1 & 1 & 1 & 1 & 277.92653 \end{array} \\ \\ \\ Mm &=& \begin{array} {rrrr} 0.0039983491 & 0.00083487726 & -0.0016554169 & -0.00076571340 \\ 0.00083487726 & 0.010062383 & 0.0027137791 & -0.0041346567 \\ -0.0016554169 & 0.0027137791 & 0.081999971 & -0.013791508 \\ -0.00076571340 & -0.0041346567 & -0.013791508 & 0.93753657 \end{array} \\ \\ \\ Gg&=&\begin{array} {rrrr} 0.0040498092 & 0.00094445419 & 0.24350752 & 0.000025246807 \\ 0.00094445419 & -0.16953880 & -0.00033124251 & 0.000014504751 \\ 0.24350752 & -0.00033124251 & -0.0010409834 & 0.0000055581782 \\ 0.000025246807 & 0.000014504751 & 0.0000055581782 & -0.000012137628 \end{array} \\ \\ \\ S &=& \begin{array} {rrrr} -5.8030040 & -3.8986593 & 14.233205 & -4.3361544 \\ -3.8986593 & -11.925784 & -1.9019096 & 1.4489718 \\ 14.233205 & -1.9019096 & 3.9412187 & -1.2246860 \\ -4.3361544 & 1.4489718 & -1.2246860 & 0.32178997 \end{array} \\ \\ \\ chk &=& \begin{array} {rrrr} -5.8030040 & -3.8986593 & 14.233205 & -4.3361544 \\ -3.8986593 & -11.925784 & -1.9019096 & 1.4489718 \\ 14.233205 & -1.9019096 & 3.9412187 & -1.2246860 \\ -4.3361544 & 1.4489718 & -1.2246860 & 0.32178997 \end{array} \\ \\ \\ err &=& \begin{array} {rrrr} 0 & . & . & . \\ . & 0 & . & . \\ . & . & 0 & . \\ . & . & . & 0 \end{array} \end{array} $
(Below is the old text of my answer, just for reference and explanation of the principle:)
I've possibly another useful reformulation, don't know whether this fits the scope of the question. I assume that A is diagonalizable and do not consider a range of convergence - it's just a sketch.
For ease of notation I introduce $ \rm\ B = transpose(A) $, small letters denote the inverse of the related capital-letters denoted matrix. $ \rm\ I $ is the identity matrix.
We have from the statement of the problem $ \rm\ S = I + B S A $
If we introduce the diagonalization of A we have $ \rm\ A = W * D * w $ where D is diagonal. Use M and m for the transpose of W and w, then the diagonalization of B is $ \rm\ B = m * D * M $
So we have $ S = I + BSA = I + mDM*S*WDw $ and thus $ MSW = MW + D MSW D $
Let's denote $ \rm\ MSW = U$ and $ \rm\ MW = t$ Then this also means: $ U = t + D U D = t + D t D + D^2tD^2 + D^3tD^3+... $ Note that t is symmetric and thus U is symmetric and $ \rm\ S = mUw $ is then symmetric as well.
Because D is diagonal this can now be described more easily for each matrix-element individually. We have, using r and c for row and column-inbdex respectively $ U_{r,c} = t_{r,c}*(1 + D_rD_c + (D_rD_c)^2 + ...) = t_{r,c}/(1-D_rD_c) $ Having the matrix U by this, S is $S = m * U * w $ [update note: I edited the first version because I'd taken A and B in the wrong order]
[edit 2]: the denominators in the description of elements of U seem to describe the range of applicability. Since geometric series with quotient q is defined also for the divergent case (except for q=1) we might discard the condition $ \varrho(A) \lt 1 $ and replace it with something like : A cannot have at the same time a value x and its reciprocal as eigenvalues