Positive moment Gram matrices produce orthonormal polynomial families

Let ${\mathcal{P}}_{p_j}$ be the real vector space spanned by the formal monomials $1,y,\dots ,y^{p_j}$. For $q,r\in {\mathcal{P}}_{p_j}$ define $⟨q,r{⟩}_{{\rho }_j}:=\int q(y)r(y)d{\rho }_j(y)$. Every product $qr$ has degree at most $2p_j$ and is integrable by [F2], so the pairing is finite; linearity of the integral makes it bilinear and symmetric. If $q(y)={\sum }_{r=0}^{p_j}{c}_r{y}^r$, then $⟨q,q{⟩}_{{\rho }_j}=c^T{G}_{j}c$, which is positive for each nonzero formal polynomial by [F1]. Thus this is an inner product on ${\mathcal{P}}_{p_j}$.

Every $q\in {\mathcal{P}}_{p_j}$ has finite squared integral by [F2], so it defines an $L^2({\rho }_j)$ class. If that class is zero, then $q=0$ almost everywhere and $\int q^{2}d{\rho }_j=0$; by [F1] its coefficient vector must be zero. Thus the class map is injective and the monomial classes are linearly independent. The formal monomials are also independent by uniqueness of their coefficient lists, so [F3] applies to this list in ${\mathcal{P}}_{p_j}$. Its residual at stage $k$ is $y^k$ minus a polynomial of degree at most $k-1$, so it has degree $k$ and leading coefficient $1$.

Each residual is nonzero in ${\mathcal{P}}_{p_j}$ by its degree and leading coefficient, so [F4] makes its norm strictly positive. The normalized residual is orthogonal to the preceding polynomials and has positive leading coefficient. At $k=0$, the residual is $1$ and its norm is $\sqrt{\int 1d{\rho }_j}=1$ because ${\rho }_j$ is a probability law; hence ${\varphi }_{j,0}=1$. This proves the stated finite orthonormal family.