Preserving explained variance while reducing dimensionality
$f$ is not defined for $\mathbb{R}^M$, so as well as reducing $x \in \mathbb{R}^N$ to $y \in \mathbb{R}^M$ you are creating a function approximator $g(y) \approx f(x)$.

I suggest using a neural network autoencoder with N dimension inputs and M dimension "bottleneck" layer. You will need to scale inputs and outputs. Using a standard auto-encoder, you will have to measure the fraction of preserved variance in $f(x)$ afterwards. The auto-encoder is not directly going to do that for you, instead it will, like PCA, attempt to use y to encode x (unlike PCA, it may do so nonlinearly). 

You could take this further if $f(x)$ is differentiable. Instead of usual mean-squared error as the loss function, you could use $\frac{1}{k}\sum_{i=1}^{k} (f(x_{i}) - f(\hat{x_i}))^2$ - that will encourage the NN to preserve variance in $x$ that matters to $f(x)$. You will have to figure out the gradient function for this analytically. Also bear in mind I have not tried this, just attempted to match your requirement to some theory.

In addition, you will now have $g(y)$ as it will be $f(\hat{x})$. You can generate any $y$ from $x$ by running the first half of the autoencoder, and can generate any $\hat{x}$ from $y$ by running second half of the autoencoder.



You might also be able to adapt t-SNE by using your function $f(x)$ to generate distances that need to be preserved when reducing dimensions.