Why does modelling t-SNE 2D coordinates with a neural network fail to produce a perfect fit? In short, I want to train a neural network to project my high dimensional objects to a 2D plane coordinates which was created using t-SNE (Rtsne package for R).
I processed 0.5 million of 128-dimensional points to this picture: 

It is just OK to me. I see the clusters. Now I train the NN with the same inputs and two outputs in an attempt to get a machine capable of mapping future samples to this plane. However, I could not get a close approximation with even very complex NN (3 layers, 1014, 512, 2 neurons), as well as with a weaker model.

So, I wonder if this is technically not possible to get exact approximation due to some randomness in the process of the t-SNE convergence, in other words, if there are contradictory examples?
Update 1 - Experiments with NN architecture:
I Z-normalized inputs to unit standard deviation and zero mean.
I normalized outputs to the range [0;1].
I use Keras for R:
Per amoeba's idea I decided to kick off with the following architecture:
# set NN

input <- layer_input(
     shape = list(as.integer(embedding_dim))
     , name = 'input'
)

common_dense1 <- layer_dense(units = 2048, activation = "elu")

x_dense2 <- layer_dense(units = 1024, activation = "elu")

x_dense3 <- layer_dense(units = 1024, activation = "elu")

x_dense4 <- layer_dense(units = 1024, activation = "elu")

x_dense5 <- layer_dense(units = 256, activation = "elu")

x_dense6 <- layer_dense(units = 1, activation = "linear")

y_dense2 <- layer_dense(units = 1024, activation = "elu")

y_dense3 <- layer_dense(units = 1024, activation = "elu")

y_dense4 <- layer_dense(units = 1024, activation = "elu")

y_dense5 <- layer_dense(units = 256, activation = "elu")

y_dense6 <- layer_dense(units = 1, activation = "linear")

tsne_model_output_x <- 
  input %>%
  #common_dense1 %>%
  x_dense2 %>%
  x_dense3 %>%
  x_dense4 %>%
  x_dense5 %>%
  x_dense6

tsne_model_output_y <- 
  input %>%
  #common_dense1 %>%
  y_dense2 %>%
  y_dense3 %>%
  y_dense4 %>%
  y_dense5 %>%
  y_dense6

tsne_output <- layer_concatenate(inputs = c(tsne_model_output_x, tsne_model_output_y))

tsne_model <- keras_model(
  inputs = input
  , outputs = tsne_output
)

keras::compile(
     tsne_model
     , optimizer = 'adam'
     , loss = 'mse'
)


## check correctness

xtrain <- array_reshape(rnorm(100 * embedding_dim), dim = c(100, embedding_dim), order = 'C')

ytrain <- array_reshape(rnorm(100 * 2), dim = c(100, 2), order = 'C')

model_train <- fit(
     tsne_model
     , xtrain
     , ytrain
     , epochs = 1
     , batch_size = 100
     , verbose = 1
)

predict(
     tsne_model
     , xtrain
)

Please note it is a parallel architecture, where each coordinate is modeled using its own layers.
I roughly estimated this model comtains 2.2 Billion of parameters.
Learning rate is 1e-3.

Final MSE: 0.013. Which means average absolute error is 0.114 (or 11% given the scale of the outputs). This is worse than I've seen before.
Next I will decrease the leraning rate to 5e-4 for another 500 epochs.
Stage 2 - smaller learning rate, larger batch size:

It seems we are stuck.
 A: I finally managed to train a neural network with desired quality (almost 100% overfitting):
# set NN

input <- layer_input(
     shape = list(as.integer(embedding_dim))
     , name = 'input'
)

dense1 <- layer_dense(units = 1024, activation = "tanh")

dense2 <- layer_dense(units = 1024, activation = "tanh")

dense3 <- layer_dense(units = 1024, activation = "tanh")

dense4 <- layer_dense(units = 2, activation = "sigmoid")

tsne_output <- 
  input %>%
  dense1 %>%
  dense2 %>%
  dense3 %>%
  dense4

tsne_model <- keras_model(
  inputs = input
  , outputs = tsne_output
)


keras::compile(
     tsne_model
     , optimizer = 'adam'
     , loss = 'mse'
)

Inputs were Z-normalized, and outputs were scaled to [0;1].
Training took about 10,000 epochs with batch size (500,000 / 64), but I varied this number during training, and learning rate that worked well was 1e-4 to 5e-5.
Generated image:

Correlation test:

Final epochs:

I could go on, but this quality is just enough.
Some thoughts about what was important:
I failed using elu activations in surpsiging way: produced image contained intersecting lines representing densed points at different angles, which made me think it was local minima and it looked weired and more like modern art. 
When I used parallel layers for x and y coordinates the image contained two orthogonal lines with dense points (like a sniper crosshair), which was also a kinf of local minima. 
I also experimented with input and output preprocessing, and output layer activation. Surprisingly, when an output layer was linearly activated, the image became quadratic.
Besides, a batch size was really important. Again, surprisingly, the larger batch size I set (let's say 500,000 / 8) the worst the converge I got. What worked well is smaller batch size, and this is really solid as I observed this numerous times.
After 3 weeks of on-and-off experiments I found that classic tanh and sigmoid and with classic preprocessing works well.
I cannot say analytically why this configuration behaves well, it is very black box, but somehow I found it intuitively.
But the funniest thing is:
Since I overfitted heavily, and I had found that overfitting starts pretty early, meaning that tSNE produces stochastic and noisy patterns to learn, the neural network prediction that I got on the new data (which are actually news texts) looks like a muddy spot as if there are no clusters. It is laughing at my efforts. Don't overfit. ) But if I did not overfit, I would get a blurry picture as well, hmm.

