I will summarize everything, for future reference.
All Lyapunov equations $AX+XA^T=B$ have a unique, symmetric solution $X=X^T$, unless there is a $\lambda\in\mathbb{C}$ such that $\lambda$ and $-\lambda$ are both eigenvalues of $A$. This holds, for instance, when $A$ is Hurwitz stable, which is a common case.
As far as I know, reducing to $B=I$ (or to $Y+Y^T=I$ with a nontrivial constraint on $Y$) does not help in its numerical solution.
There are integral representations and reformulations as a system of $n^2$ linear equations in $n^2$ unknowns using Kronecker products, but they are also worthless from the point of view of numerical solution.
Up to $n\approx 1000$, the standard algorithm (also implemented in MATLAB's
lyap) is the $O(n^3)$ Bartels-Stewart algorithm. The Hessenberg-Schur algorithm is a variation that helps for the more general Sylvester equation, but not in this case. For a general Sylvester equation $AX+XC=B$, the HS method requires one Schur decomposition of $A$ instead of one of $A$ and one of $C$, but if $C=A^T$ then the second decomposition comes for free.For values of $n$ larger than $1000$, and in general for the case of $A$ large and sparse, the standard algorithms are either rational Krylov subspace methods (MATLAB code and papers on V. Simoncini's page, http://www.dm.unibo.it/~simoncin/software.html) or the ADI method (Matlab and C code on the page of Peter Benner's group, http://www.mpi-magdeburg.mpg.de/csc). Both can handle problems where $B$ has low rank.
The case in which $A$ is large and sparse and $B$ has full rank (for instance $B=I$) is more difficult; solving it efficiently is still an active research problem. This preprint seems to go in the right direction http://www.math.cts.nthu.edu.tw/download.php?filename=683_a93bf34d.pdf&dir=publish&title=prep2012-05-005, but the methods is much newer and I have never tried it.