Just follow the Gram-Schmidt algorithm outlined there, where your inner product is the weighted inner product $$\langle u,v\rangle_w = \int_a^b u(x)v(x)w(x)\,dx.$$ Here, $w(x)$ is called the weight function.
In your context, it sounds like you want to generate the first four Laguerre polynomials, $\{L_0(x),L_1(x),L_2(x),L_3(x)\}$, by applying the Gram-Schmidt algorithm to the standard monomials $\{1,x,x^2,x^3\}$. The resulting Laguerre polynomials will form an orthogonal (or orthonormal if you include the normalization step in the Gram-Schmidt algorithm) family on $0<x<\infty$ with respect to the weight function $w(x)=e^{-x}$.
So, following the algorithm linked above (including the normalization) and using the weighted inner product above, you get (using the notation in the link):
\begin{align}
u_0(x)&=1,\\
u_1(x)&=x-{\langle x,1\rangle_w\over \langle 1,1\rangle_w}\cdot 1=x-{\int_0^\infty x\cdot 1\,e^{-x}\,dx\over \int_0^\infty 1^2\,e^{-x}\,dx}=x-1,\\
u_2(x)&=x^2-{\langle x^2,1\rangle_w\over \langle 1,1\rangle_w}\cdot 1-{\langle x^2,x-1\rangle_w\over \langle x-1,x-1\rangle_w}\cdot (x-1)\\
&=x^2-{\int_0^\infty x^2\cdot 1\,e^{-x}\,dx\over \int_0^\infty 1^2\,e^{-x}\,dx}\cdot 1-{\int_0^\infty x^2\cdot (x-1)\,e^{-x}\,dx\over \int_0^\infty (x-1)^2\,e^{-x}\,dx}\cdot (x-1)\\
&=x^2-4x+2,\\
u_3(x)&=x^3-{\langle x^3,1\rangle_w\over \langle 1,1\rangle_w}\cdot 1-{\langle x^3,x-1\rangle_w\over \langle x-1,x-1\rangle_w}\cdot (x-1)-{\langle x^3,u_2\rangle_w\over \langle u_2,u_2\rangle_w}\cdot (x^2-4x+2)\\
&=x^3-9x^2+18x-6.
\end{align}
To get the $L_k$, you just need to divide each $u_k$ above by its (weighted) norm:
\begin{align}
L_0(x)&={u_0\over \|u_0\|_w}={u_0\over \langle u_0,u_0\rangle_w^{1/2}}=1,\\
L_1(x)&={u_1\over \|u_1\|_w}={u_1\over \langle u_1,u_1\rangle_w^{1/2}}=x-1,\\
L_2(x)&={u_2\over \|u_2\|_w}={u_2\over \langle u_2,u_2\rangle_w^{1/2}}={x^2\over 2}-2x+1,\\
L_3(x)&={u_3\over \|u_3\|_w}={u_3\over \langle u_3,u_3\rangle_w^{1/2}}={x^3\over 6}-{3\over 2}x^2+3x-1.
\end{align}
and you get the first four Laguerre polynomials:
$$
\left\{1,x-1,\frac{x^2}{2}-2 x+1,\frac{x^3}{6}-\frac{3 x^2}{2}+3 x-1\right\}.
$$
Since these form an orthonormal family, you can expand $f(x)=e^{-2x}$ as follows:
$$
f(x)=\sum_{n=0}^k c_n L_n(x), \qquad c_n=\langle f,L_n\rangle_w=\int_0^\infty f(x)L_n(x)w(x)\,dx,
$$
where $k=0,1,2,3$ indicates (as you call it) the "degree" of the approximation, i.e., an expansion in terms of the first $k$ Laguerre polynomials.
For example,
\begin{align}
c_0&=\langle e^{-2x},1\rangle_w=\int_0^\infty e^{-2x}\cdot 1\cdot e^{-x}\,dx=1/3,\\
c_1&=\langle e^{-2x},x-1\rangle_w=\int_0^\infty e^{-2x}\cdot (x-1)\cdot e^{-x}\,dx=-2/9,\\
c_2&=\langle e^{-2x},\frac{x^2}{2}-2 x+1\rangle_w=\int_0^\infty e^{-2x}\cdot \left(\frac{x^2}{2}-2 x+1\right)\cdot e^{-x}\,dx=4/27,\\
c_3&=\langle e^{-2x},\frac{x^3}{6}-\frac{3 x^2}{2}+3 x-1\rangle_w=\int_0^\infty e^{-2x}\cdot \left(\frac{x^3}{6}-\frac{3 x^2}{2}+3 x-1\right)\cdot e^{-x}\,dx=-8/81.
\end{align}
Here are the $k=0,\dots,3$ approximations of $f(x)=e^{-2x}$:
We apply the Gram Schmidt orthogonalization procedure to
$$ P_1(x)=c \qquad P_2(x)=ax \qquad P_3(x)=bx^2$$
1) Take $P_1$ and normalize it
$$ ||P_1(x)||^2=\int_0^{\infty}c^2e^{-x}=c^2 $$
So $c=1 $ and $L_0=1 $
2) Take $P_2$ and subtract to it its projection on $L_0$, then normalize the result
$$P_{2|L_0} = \int_0^{\infty} (ax)e^{-x}dx =a \\P_2-P_{2|L_0}=a(x-1) \\
||P_2-P_{2|L_0}||^2 = \int_0^{\infty} a^2(x-1)^2e^{-x}=a^2$$
Now we are free to choose either $ a=1$ or $a=-1 $ , to remain consistent with your notation we choose the latter and get
$$L_1=1-x$$
3)Take $P_3$ and subtract to it its projection on $L_0$ and on $L_1$ then normalize it and get $L_2$ with very analogous procedure
Best Answer
I presume that by least squares polynomial, you mean determine coefficients $c_i$ in $$ F(x):=c_0L_0(x)+c_1 L_1(x)+c_2L_2(x)+c_3 L_3(x) $$ such that $$\|f-F\|_w^2:=\int_0^\infty |f(x)-F(x)|^2 w(x)\,dx$$ is minimized, where $w(x)=e^{-x}$.
In that case, the orthogonality of the $L_i(x)$ (and you would need to add $L_0(x)$) together with the Best Approximation Theorem say that the coefficients should be chosen according to $$c_i={\langle f,L_i\rangle_w\over \langle L_i,L_i\rangle_w}={\int_0^\infty f(x)L_i(x)e^{-x}\,dx\over \int_0^\infty L_i^2(x)e^{-x}\,dx}, \quad i=0,1,2,3. $$