An example, with some number crunching, might be intuitive and also help debunk the idea of using the initial gradient instead of $0$.
Consider the 1D problem $f(x)=x$, where $f'(x)=1$. $\beta_1=0.9$ and $\beta_2=0.999$ as usual.
The first few values of $m_t$ and $v_t$ (rounded to 4 places) are given below.
\begin{array}{c|c|c|c}
t&m_t&v_t&m_t/\sqrt{v_t}\\\hline
0&0&0&\mathrm{N/A}\\
1&0.1000&0.001000&3.162\\
2&0.1900&0.001999&4.250\\
3&0.2710&0.002997&4.950\\
4&0.3439&0.003994&5.442
\end{array}
At $t=12$, we reach a high of $m_t/\sqrt{v_t}=6.568$, and from there it descends to $1$, the "correct" value of $m_t/\sqrt{v_t}$. In other words, with these parameters we may reach step sizes roughly $6.5$ times larger than what they should be, which may be undesirable.
We can also see that initially $m_t$ and $v_t$ are very close to $0$. As dontloo shows, $m_t$ and $v_t$ are always going to start out close to the initially used value.
As Mark L. Stone comments,
What I don't understand is why the gradient at initial point is not used for the initial values of these things, and then the first parameter update. Then there would be no contamination by the initial zero values, which has to be undone. So there'd be no need for the bias correction.
Consider, however, the context in which momentum estimates are often used: stochastic/mini-batch gradient descent. It should be expected that the initial (stochastic) gradient is not an accurate estimate of the true gradient. If we truly want an accurate estimate of the gradient, then we need to have nearly equal contributions from the first few gradients.
Note then the expanded expressions for their choice of $m_t$, using $m_0=0$.
\begin{align}
m_1&=0.1g_1\\
m_2&=0.1g_2+0.09g_1\\
m_3&=0.1g_3+0.09g_2+0.081g_1\\
m_4&=0.1g_4+0.09g_3+0.081g_2+0.0729g_1
\end{align}
It is apparent that $m_t$ shares nearly the same amount of the previous several $g_t$.
Now consider setting $m_1=g_1$.
\begin{align}
m_1&=g_1\\
m_2&=0.1g_2+0.9g_1\\
m_3&=0.1g_3+0.09g_2+0.81g_1\\
m_4&=0.1g_4+0.09g_3+0.081g_2+0.729g_1
\end{align}
As expected, $g_1$ now has a ten-fold influence on $m_t$.
One could make the argument that the influence of $g_1$ in $m_t$ is rapidly diminishing, and hence largely irrelevant.
But what about $v_t$? Based on our previous example, we should expect a thousand-fold influence. Let's compare $v_t$ and $\bar v_t$, where $v_0=0$ and $\bar v_1=g_1^2$. Doing the math, here's the $\%$ influence of $g_1$ on the $v$'s for several $t$.
\begin{array}{c|c|c}
t&\%\text{ of $g_1$ in $v_t$}&\%\text{ of $g_1$ in $\bar v_t$}\\\hline
1&100\%&100\%\\
10&9.96\%&99.1\%\\
100&0.951\%&90.5\%\\
200&0.452\%&81.9\%\\
300&0.286\%&74.1\%\\
400&0.203\%&67.1\%\\
500&0.154\%&60.7\%\\
600&0.122\%&54.9\%\\
1000&0.0582\%&36.8\%\\
2000&0.0156\%&13.5\%
\end{array}
In my humble opinion, this is atrociously bad. Would I rather risk an initially inaccurate estimate of the gradient persisting in my momentum so significantly just to avoid a division by $1-\beta_2^t$? Absolutely not.
For the not-so-mathematically inclined, how does their bias correction solve all of these issues? Let's go through it, one-by-one.
At $t=12$, we reach a high of $m_t/\sqrt{v_t}=6.568$, and from there it descends to $1$, the "correct" value of $m_t/\sqrt{v_t}$. In other words, with these parameters we may reach step sizes roughly $6.5$ times larger than what they should be, which may be undesirable.
The bias correction solves this issue by rescaling $m_t$ and $v_t$ to have roughly the same magnitude as $g_t$ and $g_t^2$. How exactly? It divides the total sum by the sum of the weights of each $g_t$.
\begin{align}
\hat m_1&=\frac{0.1g_1}{0.1}\\~\\
\hat m_2&=\frac{0.1g_2+0.09g_1}{0.1+0.09}\\~\\
\hat m_3&=\frac{0.1g_3+0.09g_2+0.081g_1}{0.1+0.09+0.081}\\~\\
\hat m_4&=\frac{0.1g_4+0.09g_3+0.081g_2+0.0729g_1}{0.1+0.09+0.081+0.0729}
\end{align}
It turns out this denominator can more simply be written as $1-\beta_1^t$.
It can also be seen from the last table that by initializing $v_0=0$, we get a much more accurate momentum than initializing $v_1=g_1^2$. Indeed using $v_1=g_1^2$ actually introduces another, perhaps concerning, problem.
Furthermore, we observe that using the "initial gradient" approach is subject significantly more to the choice of $\beta$. Although $\beta=0.9$ is not so bad, when $\beta=0.999$, it can cause the initial value to persist much longer.
What should you understand from all this, intuitively? My take is that by initializing momentum to the first provided value, you become biased towards the initial value rather than biased towards $0$. In contrast, biased towards $0$ is remarkably simple, and much more intuitive, to fix.
An alternative, equivalent, formula for computing $\hat m_t$, is actually presented by dontloo, but it has some semantic drawbacks. For large $t$, we can see that $1-\beta^t\approx1$, leaving us with $m_t\approx\hat m_t$, whilst with theirs they obtain $m_t\approx\hat m_t/(1-\beta_1)$. This causes $m_t$ to be influenced by the choice of $\beta_1$ as well as lose its meaning as the momentum approximation of $g_t$. Since it's influenced by $\beta_1$, it can no longer be directly compared to $g_t$.
Though one could argue you shouldn't worry about the existence of $m_t$ and instead focus on $\hat m_t$, which is the same in both formulations, I would argue that letting $m_t$ be an approximation of $g_t$ is much more of an intuitive buildup than the latter.
For some intuition on the momentum formulas, note also the similarity between the following:
$$m_t=m_{t-1}+(1-\beta_1)(g_t-m_{t-1})$$
$$a_t=a_{t-1}+\frac1t(g_t-a_{t-1})$$
As it turns out, $a_t$ is the accumulating formula for the actual average of $g_t$. $m_t$ is then an approximation of this, where the newest gradient weighs in slightly more than the previous.