However I when thinking graphically I think that it is possible to approximate these nonlinear functions using lots of linear ( scaled and shifted) linear lines and I do not understand why this is incorrect.
It's not generally incorrect, but what you are describing there is not a single function but a piecewise function. For example:
$
f(x) = \left\{
\begin{array}{ll}
0.1 x & \quad 0 \leq x < 0.1 \\
0.2 x & \quad 0.1 \leq x < 0.2\\
0.3 x & \quad ...
\end{array}
\right.
$
This composition is non-linear, because the input space gets transformed differently depending on $x$ (because each subfunction in the example above has a different slope). There are methods to optimize such a function as well, e.g. piecewise linear regression.
An MLP without non-linearities is not a piecewise function, but a function composition:
$(f \circ g)(x) = f(g(x))$,
where $g$ and $f$ are two consecutive dense layers.
The difference is that if both functions $g$ and $f$ are linear, then $f \circ g$ is linear. This holds for any configuration of model parameters. Therefore, the model transforms any given input space $X$ into some output space $W$, where straight lines in the input space remain straight in the output space and can thus never model something non-linear.
In case, if it is possible to approximate any function using Relu activation function ( which is linear on first quadrant) , why it is not possible to approximate with a function completely linear ?
As soon as you introduce any non-linearity, the system becomes much more expressive, because each neuron now models a non-linear function. Combining the neurons (as done by a subsequent layer) can learn combinations of these non-linear functions.
And this is where the two worlds collide: If you take a ReLU activation, it allows the model to actually learn something like a piecewise linear function, because of the combinations of several linear functions with different slopes that are $< 0$ only for certain input ranges. Here is a simple example of that:
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
tf.random.set_seed(4)
np.random.seed(4)
# create a dataset
xs = np.linspace(-2, 2, 1000)
ys = xs**2 + 0.2
# create a model
m = tf.keras.Sequential([
tf.keras.layers.Dense(4, activation='relu'),
tf.keras.layers.Dense(1),
])
# train the model
m.compile(optimizer='adam', loss=tf.keras.losses.MAE)
m.fit(xs[:,np.newaxis], ys[:,np.newaxis], batch_size=8, epochs=20, verbose=None)
# plot the results
fig, ax = plt.subplots(1, 2, figsize=(15, 6))
ax[0].plot(xs, xs, label='model inputs')
ax[0].plot(xs, ys, label='ground truth')
ax[0].plot(xs, m(xs)[:,0], label='prediction')
ax[1].plot(xs, ys)
for n in range(m.layers[0].units):
neuron_act = m.layers[0](xs[:,np.newaxis])[:,n]
ax[1].plot(xs, neuron_act, color='g', label='piecewise function %d' % n)
ax[0].legend()
ax[1].legend()
ax[0].set_ylim((-1, 4))
ax[1].set_ylim((-1, 4))
The left image shows what the model has learned (green line), the right side shows the functions of the individual neurons. The more neurons you add to that first relu-layer, the better the approximation of the ground truth will be.
