You're missing a couple of quite important concepts:

• Universal approximation theorem: with enough parameters a neural network can approximate any function.
• Basically every loss function is non convex. There is this little problem in machine learning call local minima about which we like to complain a lot.

But no need to trust me, just run a simple experiment and try yourself to approximate a non convex function, like $sin(x)$ with relu:

from sklearn.neural_network import MLPRegressor
import numpy as np
import matplotlib.pyplot as plt

f = lambda x: [[x_] for x_ in x]
noise_level = 0.1
X_train_ = np.arange(0, 10, 0.2)
real_sin = np.sin(X_train_)
y_train = real_sin + np.random.normal(0, noise_level, len(X_train_))
nodes = 1000
layers = 4
regr = MLPRegressor(hidden_layer_sizes=tuple([nodes] * 4), activation="relu").fit(f(X_train_), y_train)
predicted_sin = regr.predict(f(X_train_))

plt.plot(X_train_, real_sin, label="sin target")
plt.plot(X_train_, predicted_sin, label="sin predicted")
plt.legend()
plt.show()

You'll see it's not a hard task to learn at all:

And even with a single layer (but many more hidden nodes, specifically 50k) you can get a not so good, but surely non convex solution:

You're missing a couple of quite important concepts:

• Universal approximation theorem: with enough parameters a neural network can approximate any function.
• Basically every loss function is non convex. (There is this little problem in machine learning call local minima about which we like to complain a lot :) )

But no need to trust me, just run a simple experiment and try yourself to approximate a non convex function, like $sin(x)$ with relu:

from sklearn.neural_network import MLPRegressor
import numpy as np
import matplotlib.pyplot as plt

f = lambda x: [[x_] for x_ in x]
noise_level = 0.1
X_train_ = np.arange(0, 10, 0.2)
real_sin = np.sin(X_train_)
y_train = real_sin + np.random.normal(0, noise_level, len(X_train_))
nodes = 1000
layers = 10
regr = MLPRegressor(hidden_layer_sizes=tuple([nodes] * 4), activation="relu").fit(f(X_train_), y_train)
predicted_sin = regr.predict(f(X_train_))

plt.plot(X_train_, real_sin, label="sin target")
plt.plot(X_train_, predicted_sin, label="sin predicted")
plt.legend()
plt.show()

You'll see it's not a task too hard to learn:

PS: of course this is just of a toy example, and if you decrease layers and amount of hidden units the results will become crap, but it still proves that activation surely affects, but not constrain the non linearity of the final function learned by a neural network.