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}
.
$$

2more comments