I will summarize everything, for future reference.

1. 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.

2. 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.

3. 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.

4. 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.

5. 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.

6. 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.

I will summarize everything, for future reference.

1. 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.

2. 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.

3. 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.

4. 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.

5. 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.

6. 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. EDIT: I should rather recommend this approach by a collaborator of mine https://arxiv.org/abs/1711.05493 ; it should solve exactly this problem for the case of large matrices.

I will summarize everything, for future reference.

1. 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.

2. As far as I know, reducing to $B=I$ (or to $A=I$ with some inversions of $A$) does not help in its numerical solution.

3. 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.

4. 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.

5. 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.

6. The case in which $B=I$ is more difficult, and it is 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 I the methods is much newer and I have never tried it.

I will summarize everything, for future reference.

1. 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.

2. 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.

3. 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.

4. 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.

5. 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.

6. 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.