Vtyjkyuk

$$f(x_1,x_2,...x_n):mathbbR^n rightarrow mathbbR$$
The definition of the gradient is
$$ fracpartial fpartial x_1hate_1 + ... +fracpartial fpartial x_nhate_n$$

which is a vector.

Reading this definition makes me consider that each component of the gradient corresponds to the rate of change with respect to my objective function if I go along with the direction $hate_i$.

But I can't see why this vector (defined by the definition of the gradient) has anything to do with the steepest descent.

Why do I get maximal value again if I move along with the direction of gradient?

Each component of the gradient tells you how fast your function is changing with respect to the standard basis. It's not too far-fetched then to wonder, how fast the function might be changing with respect to some arbitrary direction? Letting $vec v$ denote a unit vector, we can project along this direction in the natural way, namely via the dot product $textgrad( f(a))cdot vec v$. This is a fairly common definition of the directional derivative.

We can then ask in what direction is this quantity maximal? You'll recall that $$textgrad( f(a))cdot vec v = |textgrad( f(a))|| vec v|textcos(theta)$$

Since $vec v$ is unit, we have $|textgrad( f)|textcos(theta)$, which is maximal when $cos(theta)=1$, in particular when $vec v$ points in the same direction as $textgrad(f(a))$.

Consider a Taylor expansion of this function,
$$f(bf r+bfdelta r)=f(bf r)+(nabla f)cdotbfdelta r+ldots$$
The linear correction term $(nabla f)cdotbfdelta r$ is maximized when $bfdelta r$ is in the direction of $nabla f$.

Other answers are correct in using the directional derivative to show that the gradient is the direction of steepest ascent/descent. However, I think it is instructive to look at the definition of the directional derivative from first principles to understand why this is so (it is not arbitrarily defined to be the dot product of the gradient and the directional vector).

Let $f(mathbfx):mathbbR^n rightarrow mathbbR$. The partial derivatives of $f$ are the rates of change along the basis vectors of $mathbfx$:

$textrmrate of change along mathbfe_i = lim_hrightarrow 0 fracf(mathbfx + hmathbfe_i)- f(mathbfx)h = fracpartial fpartial x_i$

Each partial derivative is a scalar. It is simply a rate of change.

The gradient of $f$ is then defined as the vector:

$nabla f = sum_i fracpartial fpartial x_i mathbfe_i$

We can naturally extend the concept of the rate of change along a basis vector to a (unit) vector pointing in an arbitrary direction. Let $mathbfv$ be such a vector, i.e., $mathbfv = sum_i alpha_i mathbfe_i$ where $sum_i alpha_i^2 = 1$. Then:

$textrmrate of change along mathbfv = lim_hrightarrow 0 fracf(mathbfx + hmathbfv) - f(mathbfx)h$

Again, this quantity is a scalar.

Now, it can be proven that if $f$ is differentiable at $mathbfx$, the limit above evaluates to: $(nabla f) cdot mathbfv$. This is a dot product of two vectors, which returns a scalar.

We know from linear algebra that the dot product is maximized when the two vectors point in the same direction. This means that the rate of change along an arbitrary vector $mathbfv$ is maximized when $mathbfv$ points in the same direction as the gradient. In other words, the gradient corresponds to the rate of steepest ascent/descent.

The question you're asking can be rephrased as "In which direction is the directional derivative $nabla_hatuf$ a maximum?".

Assuming differentiability, $nabla_hatuf$ can be written as:

$$nabla_hatuf = nabla f(textbfx) cdot hatu =|nabla f(textbfx)||hatu|cos theta = |nabla f(textbfx)|cos theta$$

which is a maximum when $theta =0$: when $nabla f(textbfx)$ and $hatu$ are parallel.

Each component of the derivative
$$
fracpartial fpartial x_1 ... fracpartial fpartial x_n$$
tells you how fast your function is changing with respect to the standard basis.
It's now possible to make a basetransformation to an orthogonal base with $ n-1 $ base Directions with $0$ ascent and the gradient direction. In such a base the gradient direction must be the steepest since any adding of other base directions adds length but no ascent.

For a 3 dimensional Vector space the base could look like this
$$ left( left( beginmatrix partial x_2 \ -partial x_1 \ 0 endmatrix right) left( beginmatrix partial x_1 \ partial x_2 \ -dfrac(partial x_1)²+(partial x_2)²partial x_3 endmatrix right) left( beginmatrix partial x_1 \ partial x_2 \ partial x_3 endmatrix right) right) $$ By complete induction it can now be shown that such a base is constructable for an n-Dimensional Vector space.
$$ left( left( beginmatrix partial x_2 \ -partial x_1 \ 0 \ 0 endmatrix right) left( beginmatrix colorbluepartial x_1 \ partial x_2 \ -dfrac(partial x_1)²+(partial x_2)²partial x_3 \ 0 endmatrix right) left( beginmatrix colorbluepartial x_1 \ partial x_2 \ colorgreenpartial x_3 \ -dfrac(partial x_1)²+(partial x_2)²+(partial x_3)²partial x_4 endmatrix right) left(beginmatrix colorbluepartial x_1 \ partial x_2 \ colorgreenpartial x_3 \ colororangepartial x_4 endmatrix right) right) $$ One can see here that the first Basevector demands the first 2 Elements of the following Basevectors to be $partial x_1$ & $partial x_2$ because of the orthogonal condition,
similarly the 2nd vector demands all the 3rd elements of the following vectors to be $partial x_3$
as does the 3rd vector for the 4th element them being $partial x_4$.

If another dimension is added the n+1 Element of the n$th$ Vector needs to be $$-dfrac(partial x_1)²+...+(partial x_n)²partial x_n+1$$ to meet the $0$ ascension condition which in turn forces the new n+1$th$ Vector to be of the form $$left(beginmatrixpartial x_1 \ ... \ partial x_n+1endmatrixright)$$ for it to be orthogonal to the rest.

Let $vec v$ be an arbitrary unit vector. Then the change of $f$ by moving in the direction of $v$, starting in point $a$, is given by $grad( f(a)) cdot vec v$. We want to find a $vec v$ for which this inner product is maximal.
For the inner product we have the Cauchy–Schwarz inequality $vec a cdot vec b leq |vec a||vec b|$.
Now the equality holds when $vec v = lambda ; grad(f(a))$, for some $lambda in mathbbR$.

Let $v=fracss$ be a unit vector and assume that $v$ is a descent direction, i.e. $v^Tnabla f(x) <0$. Then $f(x+lambda v)$ as a function of $lambda$, describes how this function changes along the direction $v$.

The rate of descent at $x$ along $v$ is given by:
$$ fracdd lambdaf(x+lambda v)|_lambda=0 = v^T nabla f(x) =fracs^Tsnabla f(x) equiv fracs^Tsg$$
So we want to find the maximum of this quantity as a function of $s$.
Differentiating the above wrt $s$ and setting it equal to zero, we get (noting that $nabla_s|s| =fracss$): $g=(g^T v)vequiv av$.

Taking the Euclidean norm: $|g|=|a||v|=|a| Rightarrow a=pm|g|$.

We choose the minus sign to satisfy that $v$ is descent. Hence the direction of the steepest descent is
$$ v= dfrac1ag = -dfracgg$$

Just want to further clarify why the gradient provides the steepest ascent (instead of descent) here. Any differentiable $f$ can be approximated by the linear tangent plane, i.e., $$f(mathbfx + h mathbfv) = f(mathbfx) + h , nabla f(mathbfx)^T mathbfv $$ as $h rightarrow 0$ for any unit-length direction $mathbfv$ with $parallel mathbfv parallel =1.$ As $h downarrow 0$, consider the amount of change
$$
f(mathbfx + h mathbfv) - f(mathbfx) = h , left , nabla f(mathbfx)^T mathbfv right
~~in~~ left[ - h , parallel nabla f(mathbfx) parallel, ~ h , parallel nabla f(mathbfx) parallel right]
$$
by Cauchy-Swcharz inequality, which reaches its maximum (increase) $(h , parallel nabla f(mathbfx) parallel)$ when $mathbfv = nabla f(mathbfx) / parallel nabla f(mathbfx) parallel$ and its minimum (i.e., maximum decrease) $ (-h , parallel nabla f(mathbfx) parallel) $ if $ mathbfv= - nabla f(mathbfx)/parallel nabla f(mathbfx) parallel$ (the negative gradient direction).