How do the infinitesimal generators generate the whole compact Lie group? Suppose that $G$ is a connected compact Lie group, $g(a)\in G$ is its element, and $a=(a_1,a_2,...,a_r)$ is the parameter of the group element $g(a)$. Then the infinitesimal generator of the group can be defined as 
\begin{equation}
X_i=\frac{\partial{g(a)}}{\partial{a_i}}\bigg|_{a=0} 
\end{equation}
Then an infinitesimal element $g(\delta a)=E+\sum_i \delta a_i X_i$, in which
$\delta a\in O({\bf0})$, $E=g({\bf0})$ is the identity element and $O({\bf0})$ is the neighborhood of the ${\bf 0}$ point in parameter space. We know that the parameter of the product of two infinitesimal elements is the sum of the parameters of the two infinitesimal elements (to the first order of the infinitesimal parameter), that is
$$
g(\delta a)g(\delta b)=(E+\sum_i \delta a_i X_i)(E+\sum_i \delta b_i X_i)=E+\sum_i (\delta a_i+\delta b_i) X_i=g(\delta a+\delta b). \tag{1}
$$
Then, in every book I have read, it says that, any finite element $g(a)$ can be generated by $\{X_i\}$: let $\delta a=a/N$,
\begin{align}
g(a)=g(N\frac{a}{N})=g(N\delta a)&=g(\delta a+\delta a+...+\delta a)\\
 &=*=g(\delta a)g(\delta a)...g(\delta a)=g(\delta a)^N\\
 &=(E+\sum_i\delta a_iX_i)^N=(E+\frac{\sum_i a_iX_i}{N})^N\\
 &=\exp(\sum_i a_iX_i)~~[\text{let } N\rightarrow \infty] \tag{2}
\end{align}
In the above proof, the key step is the $(=*=)$ step, which uses the relation in Eq.(1) $g(\delta a+\delta b)=g(\delta a)g(\delta b)$. But this relation is ONLY valid for infinitesimal elements and between two infinitesimal elements, it can not be used to derive the $(=*=)$ step directly because  both $a$ is finite and there are $N$ elements here.
My questions are:


*

*Since the ($=*=$) step is wrong, is the conclusion Eq.(2) also wrong? (not limited to this proof)

*If the the conclusion Eq.(2) is right, how to prove it in a right way?


I think I have a counterexample. Using the same logic as the above proof, we consider a scalar function $f(x)$ with single parameter. It has the property $f(0)=1$. Then the "infinitesimal generator" of $f(x)$ is just its derivative at 0, let it be
$$
t=\frac{\partial f(x)}{\partial x}\bigg|_{x=0}=f'(0)
$$
Then for any infinitesimal parameter $\delta x$ and $\delta y$, to the first order, we have $f(\delta x)=f(0)+f'(0)\delta x=1+\delta x\cdot t$ and $f(\delta y)=1+\delta y\cdot t$, and 
$$
  f(\delta x)f(\delta y)=(1+\delta x\cdot t)(1+\delta y\cdot t)=1+(\delta x+\delta y)t=f(\delta x+\delta y)
$$
Then, if the logic in Eq.(2) is right, we can also derive the conclusion that: "For any finite $x$, there is the relation $f(x)=\exp(tx)$". But this conclusion is obviously wrong! Because here $f(x)$ can be any function satisfying $f(0)=1$, e.g. $f(x)=1+\sin(x)$. Obviously, $1+\sin(x)\ne\exp(tx)$ !  
Only when $f(x)$ satisfies $f(x+y)=f(x)f(y)$ for any $x$ and $y$ (not limited to infinitesimal parameters), we can derive $f(x)=\exp(tx)$.
 A: This question stems from a sequence of misconceptions about Lie groups, I will address only some of these. 


*

*
Suppose that $G$ is a connected compact Lie group, $g(a)\in G$ is its element, and $a=(a_1,a_2,...,a_r)$ is the parameter of the group element $g(a)$.

This is not a meaningful sentence. A Lie group $G$ is a differentiable manifold equipped with a group structure satisfying certain properties. In particular, elements $g$ of $G$ do not come equipped with parameters $a=(a_1,...,a_n)$. One can say that for every $g\in G$ there exists a chart $(U,\phi)$ of the smooth atlas on $G$ (here $g\in U$, $\phi: V\subset R^n\to U$ is a diffeomorphism defined on an open subset of $R^n$). One can also say that for elements of $G$ sufficiently close to the identity, one can use a chart $(U,\phi)$ given by the exponential map of $G$. Attempting to do a calculation in (1) or/and (2) in your question will utterly fail for arbitrary charts. However, (1) can be made into a meaningful statement is we use  the exponential map as a chart. 


*
Then the infinitesimal generator of the group can be defined as 
  \begin{equation}
X_i=\frac{\partial{g(a)}}{\partial{a_i}}\bigg|_{a=0} 
\end{equation}

What you are saying here is that we can take a basis $X_1,...,X_n$ in the tangent space $T_eG$ and declare its elements to be infinitesimal generators of $G$. 


*
Then an infinitesimal element $g(\delta a)=E+\sum_i \delta a_i X_i$, in which
  $\delta a\in O({\bf0})$, $E=g({\bf0})$ is the identity element and $O({\bf0})$ is the neighborhood of the ${\bf 0}$ point in parameter space.

The mathematical translation of this sentence is that the vectors $X_1,...,X_n$ span the Lie algebra ${\mathfrak g}$, which is a part of the definition of a basis of a vector space. The staff about neighborhoods is irrelevant here. 


*
We know that the parameter of the product of two infinitesimal elements is the sum of the parameters of the two infinitesimal elements (to the first order of the infinitesimal parameter), that is
  $$
g(\delta a)g(\delta b)=(E+\sum_i \delta a_i X_i)(E+\sum_i \delta b_i X_i)=E+\sum_i (\delta a_i+\delta b_i) X_i=g(\delta a+\delta b). $$

To make sense of this statement you either have to assume that $G$ is realized as a subgroup of $GL(N, {\mathbb R})$ for some $N$ (this is a nontrivial theorem that every compact connected Lie group is linear) or you can say something along the lines: Equip $G$ with a Riemannian metric. Then for $X, Y\in {\mathfrak g}$ with sufficiently small norm, we have
$$
dist(\exp(X+ Y), \exp(X) \exp(Y))= O((||X||+||Y||)^2), 
$$
where $dist$ is the Riemannian distance function on $G$. 


*The next sentence



Then, in every book I have read, it says that, any finite element $g(a)$ can be generated by $\{X_i\}$: let $\delta a=a/N$,
  \begin{align}
g(a)=g(N\frac{a}{N})=g(N\delta a)&=g(\delta a+\delta a+...+\delta a)\\
 &=*=g(\delta a)g(\delta a)...g(\delta a)=g(\delta a)^N\\
 &=(E+\sum_i\delta a_iX_i)^N=(E+\frac{\sum_i a_iX_i}{N})^N\\
 &=\exp(\sum_i a_iX_i)~~[\text{let } N\rightarrow \infty] \tag{2}
\end{align}

is meaningless, so no wonder that you have hard time interpreting it and running into contradictions. As I said in Part (1), you cannot use some random chart as your parameter for (1) to make sense; you have to use the exponential map. Then (2) becomes a "proof" of what you already know, namely, $g=\exp(X)$ for some $X\in {\mathfrak g})$ (which you want for some reason to write as a linear combination $X=\sum_i a_i X_i$ of the basis elements) provided that $g$ is sufficiently close to the identity in $G$. or, it becomes a "proof" of a nontrivial mathematical result, namely that for compact connected Lie groups the exponential map is surjective, hence, for every $g\in G$ there exists $a\in {\mathbb R}^n$ such that $g=g(a)$, meaning that
$$
g= \exp(\sum_{i=1}^n a_i X_i).
$$
The computation you quoted cannot prove this fact, in particular, since it does not use compactness of $G$ and for connected noncompact (even linear) Lie groups the exponential map is not, in general, surjective.  At best, this computation "verifies" the fact that for every $X\in {\mathfrak g}$ and every $k\in {\mathbb N}$,
$$
\exp(X)= \left(\exp(\frac{1}{k}X)\right)^k, 
$$
which is true even if $G$ is noncompact. The proof given in (2) is incomplete. A complete proof would go along the following lines:
First verify that if ${\mathfrak a}\subset {\mathfrak g}$ is a commutative Lie subalgebra then for any two elements $X, Y\in {\mathfrak a}$, one has
\begin{equation}\label{*}
\exp(X+Y) =\exp(X) \exp(Y)  \tag{3}
\end{equation}
(no error term!). This can be proven for instance via Baker–Campbell–Hausdorff formula. 
Remark. The proof of this equation is easier if you assume that $G$ is (compact) linear and that $Y$ is a scalar multiple of $X$. Then (3) 
follows from the fact that $X$ is diagonalizable (over the complex numbers) and that (3) clearly holds for diagonal matrices.  
Then use this formula repeatedly for real multiples of $X$. This will work, but will shed no light on your question "how infinitesimal generators generate the Lie group". 
The best proof of surjectivity of the exponential map for compact connected Lie groups that I know uses Riemannian geometry: First equip $G$ with a biinvariant Riemannian metric, then verify that the Lie-exponential map at the identity coincides with the Riemannian exponential map. Lastly, use a corollary of Hopf-Rinow theorem that the exponential map is surjective for every compact connected Riemannian manifold.  
A: Formulas (1) and (2) are fine, and so is (*), properly interpreted. But you may have been misreading the otherwise sensible texts you dealt with; you might consider some texts recommended by the Greek chorus of comments above. 
Let's start from your false counterexample, which goes to the heart of your misconception. It is, of course, Lie's celebrated advective flow. (Beware of the perverse inverted use of variables from here!)
Any continuous one-parameter Abelian group is equivalent to the group of translations. Take your
$$
f(z+\delta x) \equiv g(\delta x) ~ f(z)= \left (1+\delta x  \frac{\partial}{\partial z}\right )f(z) +O(\delta x ^2),
$$
to leading order in $\delta x$. Here, g represents the finite, not the infinitesimal group element: it is only approximated by the binomial. Note the generator t does not involve a specific f(x), it involves all of them.
You may now translate further
$$
f(z+\delta x + \delta y)= g(\delta y)f(z+\delta x)=g(\delta x)g(\delta y)f(z)=g(\delta x +\delta y)f(z).
$$
Composing N of these group elements amounts to Taylor expanding around an ever-shifting argument, not just z, as you apparently are reading the group composition to mean. The group property (*) always works, not just for infinitesimal argument. The Nth power of the binomial is only the dominant term of the Nth power of g, corrected by $O(1/N)$.
For infinite N, you have your exponential, Lagrange's translation operator which summarizes the Taylor expansion compactly,
$$
g(x)f(z)=e^{xt} f(z)=e^{x \frac{\partial}{\partial z }} f(z)= f(z+x).
$$
You may now set z=0; I kept it to emphasize t is an operator--a gradient of the variable of all f. Thus, $g(x)f(0)=f(x)$.
This could all appear like a triviality, except, in principle, you (Lie) have now solved the problem for arbitrary coordinates and flows.
That is, if, instead of the gradient, you had a messier (reasonable) function, 
$$
t=\beta(z) \partial_z,
$$
You could change variable z to canonical coordinates h(z) s.t. 
$$
h'(z)=1/\beta(z), 
$$
so that 
$$
t=\frac{\partial}{\partial h} ,
$$
and you'd use this "Abel function" to translate your f through functional conjugacy. (Read up in WP, if interested, but no matter.)
Back to you question. Since you are referring to physics texts, let's stick to n×n matrices for the generators X. For one generator t and one parameter a (a unique sum as you have), you just repeat the above. This is, in fact, equivalent to the definition of the matrix exponential. 
For two, you utilize the Trotter formula, well, actually the whole enchilada, the CBH expansion formula, underlying Lie's third fundamental theorem, to see that individual exponentials compose to a generic exponential of Lie Algebra elements.
When in doubt, check with familiar elementary rotations on a sphere (our globe!), using Pauli or 3x3 matrices, which I'm sure you've mastered before moving on to recondite Lie Theory texts.
