Kevin Smith wrote:
I'm almost certain the answer to my question has been known for centuries, but I don't know where to find it...
You're right! It goes back to Euler. You're looking for a variation on Euler's infinite product representation algorithm (EIPRA). EIPRA takes
as input a sequence $b_n$ and outputs a sequence $a_n$ such that
$$1+\sum_{n=1}^{\infty}b_n x^n = \prod_{n=1}^{\infty}\left(1-x^n\right)^{-a_n},$$
whereas you are looking for an algorithm that takes as input a sequence $b_n$ and outputs a sequence $a_n$ such that
$$1+\sum_{n=1}^{\infty}b_n x^n = \prod_{n=1}^{\infty}\left(1-a_nx^n\right)^{-1}.$$
Fortunately, the key to both algorithms is to take the logarithmic derivative of both sides of the
equation, which transforms the infinite product into a Lambert series.
[I see from the comments
that you already considered logarithmic differentiation, so you must be looking for
something more explicit than the recurrence that follows.]
Taking the logarithmic derivative of both sides of your equation (1) (and multiplying by $x$),
$$
\sum_{n=1}^{\infty}\frac{a_n n x^n}{1-a_n x^n}=
\frac{\sum_{n=1}^{\infty}b_n n x^n}{ 1+\sum_{ n=1}^\infty b_n x^n }
:= \sum_{n=1}^\infty c_n x^n .
$$
The sequence $c_n$ is easily determined by the sequence $b_n$, and the sequence $a_n$ is then determined from the sequence $c_n$ as in EIPRA.
Expanding $\frac{1}{1-a_n x^n}$ as a geometric series,
$$
\sum_{m=1}^\infty \sum_{j=1}^\infty m a_m^j x^{m j}
= \sum_{n=1}^\infty c_n x^n ,
$$
hence, $a_n$ is defined recursively by
$$
n a_n
= c_n - \sum_{m|n,m\ne n} m a_m^{\frac n m}
.
$$