With $V$ American option value, $H$ holding (aka continuation) value, and $B$ bank account value, we have:
$$V_N(S_N) = (K-S_N)^+$$
and for $i$ backwards from $N-1$ down to 0, we have:
$$ H_i(S_i) = \mathbf{E}\left[B_{i}B_{i+1}^{-1}V_{i+1}(S_{i+1})|S_i\right]$$
$$ V_i(S_i) = \max (K-S_i, H_i(S_i)) $$
The algorithm result is $V_0(S_0)$.
The different regression functions $r_i$ (linear combinations of basis functions) give the conditional expectations needed at every step $i$:
$$\mathbf{E}\left[B_{i}B_{i+1}^{-1}V_{i+1}(S_{i+1})|S_i\right] = r_i(S_i)$$
Optimal stopping time index $\eta$ taking values in $\{1,..., N \}$ is defined as:
$$ \eta = \min \{k\geq 1 \mid K-S_k \geq r_k(S_k) \} \wedge N$$
Edit: What is the exercise boundary in the probabilistic (Monte Carlo) framework?
When solving for the optimal stopping time index $\eta$ above, one needs to introduce the probability state set $\Omega =\{\omega^1,...,\omega^J\}$, which in Monte Carlo context represents the labels of the primitive market variables fully (all the way to the end of the financial contract) simulated paths, as $\eta$, like all variables here, are functions on it:
$$ \eta (\omega^j)= \min \: \{k\geq 1 \mid K-S_k(\omega^j) \geq r_k(S_k(\omega^j)) \} \wedge N$$
The underlying price on path $\omega_j$ where one stops is
$$ S_{\eta{(\omega_j)}} (\omega_j).$$
If we introduce a new variable $\Gamma: \{1,...,n\}\times \Omega \rightarrow \{0,1\}$ defined as:
$$ \Gamma (i, \omega) = 1 {\rm \: if \:} i > \eta(\omega), $$
and $0$ otherwise, then the (topological) boundary of the set
$$\Gamma^{-1}(1) = \{(i,\omega) | i > \eta(\omega) \} $$ is called the exercise boundary.
Edit2: If you are referring to the PDE framework, then indeed the exercise boundary is defined as the deterministic curve $S^{\rm opt}_k$
$$ S^{\rm opt}_k = \inf \: \{x | (K - x)^+ = V_k(x) \} $$
Edit3: If one uses the conditional density $\phi_i(y\mid x)$ (transition from $S_i$ to $S_{i+1}$) to compute the conditional expectation for holding value by integration (rather than using regression function $r_i$), we have
$$ H_i(x) = \mathbf{E}\left[B_{i}B_{i+1}^{-1}V_{i+1}(S_{i+1})|S_i =x\right]$$
$$ = B_{i}B_{i+1}^{-1}\int_0^\infty \max (K-x, H_{i+1}(x))\phi_i(y\mid x)dx $$
The integral calculation can be improved if we split its integration domain at $x^*$, the root of the equation:
$$ K-x = H_{i+1}(x)$$
