Dimensionality reduction with least distance distortion Question: Could I find a dimensionality reduction algorithm without or with minimal distance (cosine) distortion?
Background: I would like to visualize in 2D a sample of news texts for which I also apply clustering.
I present the texts as vectors (like text2vec), but it is desirable to also project the vectors to 2 dimensions. Moreover I consider apply clustering on the 2 dimensions for speedy execution.
My initial guess is using t-SNE with careful tuning of epsilon.
Are there more robust algorithms which do not distort much the distance (cosine) in multidimensional space, thus making clear pictures?
Update.
After several tries I got the following mapping (and clustering based on that):
This is a daily snapshot of Russian news by three major news agencies.
People liked the look, but what I worry about is that in the middle is a huge all-in-one cluster where all kinds of topics have been mixed. However, I know that there are many smaller topics every day that consist of maybe 3-5 news, but the meaning of them is quite specific.
I played quite some time with the parameters of t-SNE:
## rt-sne reduction

mat_dat <- data.matrix(text_centroids[, !c('message_id'), with = F])

rtsne_obj <- Rtsne(X = mat_dat
                   , dims = 2
                   , initial_dims = 100
                   , perplexity = 30
                   , theta = 0.3
                   , check_duplicates = FALSE
                   , pca = TRUE
                   , partial_pca = FALSE
                   , max_iter = 1000
                   , verbose = T
                   , is_distance = FALSE
                   , Y_init = NULL
                   , pca_center = TRUE
                   , pca_scale = TRUE
                   , normalize = TRUE
                   , momentum = 0.5
                   , final_momentum = 0.8
                   , eta = 200
                   , exaggeration_factor = 12
                   , num_threads = 1)

text_centroids[, Dim1 := rtsne_obj$Y[,1]]
text_centroids[, Dim2 := rtsne_obj$Y[,2]]

ggplot(data = text_centroids) + 
     geom_point(aes(x = Dim1, y = Dim2), color = 'blue', alpha = 0.1, size = 4) +
     theme_minimal()


## cluster with dbscan

clustering_dat <- rtsne_obj$Y[, 1:2]

     dbscan_knn <- frNN(x = clustering_dat, eps = 2, sort = TRUE, search = "kdtree", bucketSize = 10,
                        splitRule = "suggest", approx = 0)

     dbscan_obj <- dbscan(x = dbscan_knn
                          , weights = NULL
                          , borderPoints = F
                          )

     table(dbscan_obj$cluster)

but could not get any better.
Update 2.
Measuring Pearson's correlation coefficient between the vector of Euclidean distances between full-dimensional points (right triangle of the dist matrix) and t-SNE dimensinal points yielded:
> cor.test(original_text_dist, tsne_dist)

    Pearson's product-moment correlation

data:  original_text_dist and tsne_dist
t = 1095.5, df = 2388200, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5774797 0.5791678
sample estimates:
      cor 
0.5783244 

Hm, it looks that the distances were distorted although the distrortion was not so catastrophic.
I think I could try multidimensional scaling next.
Hmm, it looks like cmdscale in R stats:: does significantly better job at preserving the distances:
> mds <- cmdscale(original_text_dist)
> mds_dist <- dist(mds, diag = F, method = "euclidean", upper = F)
> original_text_dist <- as.vector(original_text_dist)
> mds_dist <- as.vector(mds_dist)
> cor.test(original_text_dist, mds_dist)

    Pearson's product-moment correlation

data:  original_text_dist and mds_dist
t = 1834, df = 2388200, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.7641945 0.7652476
sample estimates:
      cor 
0.7647216 


 A: The one possible approach could be a relatively new trimap method, based on triplets (in contrast with usual pair-wise distances). It's explicitly optimize to keep the original form (in contrast to t-sne), while making better clusters then PCA.

https://github.com/eamid/trimap
A: Let me answer myself.
First of all, my criterion of choosing any method for dimensionality reduction is visual inspection of a 2D mapping with a layer of clusters and their description. Judging by that, I like t-SNE more given it has mapped some of my texts into clear clusters (topics).
MDS and PCA that I also tried did not give a clustering look to the mapping.
Another option that I explored is a neural network designed specifically (by myself) for this task. Take a look and comment if you think it is worth it.
First, I create input vector which is a concatenation of two objects in the k-dimensonal space.
Then I build a multilayer dense network with bottleneck. The k-dimesnional input vector is being convolved to just 2 scalars.
Notice that the layers that convolve k dimensions to 2 dimensions are shared between both input vectors.
After that I include a stack of dence layers that should approximate euclidean distance (calculated between k-dimensional objects) using just 2-dimensional objects stemming from the bottleneck part.
As a result the NN leans to approximate the euclidean distance from first calculating an input of 2 + 2 scalars denoting the preojection on a 2D plane.
## build keras model for dimension reduction -------------------------

input1 <- layer_input(shape = c(dim(dat_train)[2] / 2), name = 'input1')

input2 <- layer_input(shape = c(dim(dat_train)[2] / 2), name = 'input2')

shared_layer2 <- layer_dense(units = 128, activation = "relu")
shared_layer3 <- layer_dense(units = 32, activation = "relu")
shared_layer4 <- layer_dense(units = 2, name = 'shared_layer4')

text_1_model <- input1 %>%
     shared_layer2 %>% 
     shared_layer3 %>% 
     shared_layer4

text_2_model <- input2 %>%
     shared_layer2 %>% 
     shared_layer3 %>% 
     shared_layer4

main_output <- layer_concatenate(c(text_1_model, text_2_model)) %>%
     layer_dense(units = 32, activation = "relu") %>%
     layer_dense(units = 16, activation = "relu") %>%
     layer_dense(units = 1)

convolving_model <- keras_model(
     inputs = c(input1, input2), 
     outputs = main_output
)

summary(convolving_model)

convolving_model %>% compile(
     loss = "mse",
     optimizer = optimizer_adam(),
     metrics = list("mean_squared_error")
)

epochs <- 500
batch_size <- dim(train.x.1)[1]

# Fit the model and store training stats

history <- convolving_model %>% fit(
     x = list(train.x.1, train.x.2)
     , y = train.y
     , batch_size = batch_size
     , epochs = epochs
     , validation_split = 0.3
)

plot(history)

After the NN has converged, I take the output of a hidden layer that is the 2D plane projections of the two input vectors. Use them as your dimensionality reduction matrix of coordinates.
## get hidden representations

layer_name <- 'shared_layer4'

intermediate_layer_model1 <- keras_model(
                              inputs = convolving_model$input
                              , outputs = get_output_at(shared_layer4, 1)
                                        )

intermediate_layer_model2 <- keras_model(
     inputs = convolving_model$input
     , outputs = get_output_at(shared_layer4, 2)
)

intermediate_output1 <- as.data.table(predict(intermediate_layer_model1, list(train.x.1, train.x.2)))

intermediate_output2 <- as.data.table(predict(intermediate_layer_model2, list(train.x.1, train.x.2)))

calc_coords <- cbind(intermediate_output1, intermediate_output2)

Surprisingly, the correlation of the euclidean distances between the calculated 2D coordinates and the original k-dimensional distances is a whopping 0.85 (if the NN converged well).
colnames(calc_coords) <- c(
     'x1'
     , 'y1'
     , 'x2'
     , 'y2'
)

calc_coords[, calc_euclid := sqrt((x1 - x2)^2 + (y1 - y2)^2)]

cor.test(train.y, calc_coords[, calc_euclid])

