For multi-valued functions-such as when the conditional distribution on targets given the inputs is multi-modal-standard regression approaches are not always desirable because they provide the conditional mean. Modal regression approaches aim to instead find the conditional mode, but are restricted to nonparametric approaches. Such approaches can be difficult to scale, and make it difficult to benefit from parametric function approximation, like neural networks, which can learn complex relationships between inputs and targets. In this work, we propose a parametric modal regression algorithm, by using the implicit function theorem to develop an objective for learning a joint parameterized function over inputs and targets. We empirically demonstrate on several synthetic problems that our method (i) can learn multi-valued functions and produce the conditional modes, (ii) scales well to high-dimensional inputs and (iii) is even more effective for certain uni-modal problems, particularly for high frequency data where the joint function over inputs and targets can better capture the complex relationship between them. We conclude by showing that our method provides small improvements on two regression datasets that have asymmetric distributions over the targets.
This paper considers the regression problem in scenarios in which the conditional distribution of the response variable y given the input x is multimodal. For artificially constructed datasets, in which the the conditional modes are known, performance is assessed by the RMS distance to the mode closest to each prediction.  For real-world datasets, performance is in terms of RMSE and MAE.
The authors propose to learn a function f(x,y) that takes both an input x and a proposed output y and produces a value of 0 when the output y is correct for the input (i.e a conditional mode of p(y|x)), and hopefully a nonzero output when y is incorrect.  To predict on a point x, one searches for a y that makes f(x,y)=0.  For this to work, we will need f(x,y)≠0 for incorrect y. To achieve this, the authors take inspiration from the implicit function theorem and additionally seek an f with ∂f/∂y=-1 at the modes for each x.  For implementation, they seeks an f that minimizes the objective f²(xᵢ,yᵢ) + ( ∂f(xᵢ,yᵢ)/∂y - 1 )² over the training data.  To predict y for a given x, they find argmin_y f²(x,y) + ( ∂f(x,y)/∂y - 1 )² by grid search over the range of y, which they assume to be known. 
I think this is a very interesting and novel approach to regression.  It seems to work well in artificial examples and competitively on some real world examples.  
My main concern is that, although the objective function is fairly intuitive, it is not  clear to me what the properties are of the population minimizer (i.e. in the infinite data case).  Equation (1) defines modes as local maxima of the density, and for evaluation purposes on the artificial data, we look at deviation from the closest conditional mode.  With infinite data, do we expect the f to converge to something that's 0 at all the conditional modes, or just at the ones with the largest density?  What do we even want to happen?  Can we use this approach to predict all the modes for each x?  While I consider all of these interesting questions, after some consideration, I don't think that missing answers to these questions should delay publication.  
I recommend accepting this paper, assuming reasonable responses can be provided to the questions and issues below.  It's a novel and promising approach to an interesting problem in regression.  It leaves some open questions for further theoretical analysis, but that seems ok to me.
Specific questions:
- In section 3, second sentence, I think you want to write either the minimum of f^2 or a zero of f?
- Page 3, not sure what is meant by "should be described by the same mapping g_j(x)."
- In the implicit function theorem, it ends with a capital F, but I think you want f.
- When you say "this condition is only local and may not encourage a minimal set of conditional modes" -- I guess you mean a minimum set of zeros? Or predicted conditional modes?
- In Equation 5, you write a full l_2 norm for what I believe is just a scalar -- how about parenthesis instead?
- When you make your y prediction using the argmin... it seems like you could also minimize without the partial derivative term, as in your original explanation of f.  Did you compare the performance with and without the partial derivative term?
- In comparing to the MDN network, it's interesting that MDN needs more mixture components than the true distribution, but that doesn't really seem like a point of criticism.  Neural networks themselves seem to work better with far more parameters than we think they should need.  To be fair, I'd suggest hyperparameter searching over a broader range of mixture components for MDNs, maybe up to 10.  
- In your section "verify the learned error distribution", you say that "it is obvious that the error function indeed looks like a Gaussian distribution and the one trained without noise 4(a) shows an extremely small variance."  Two comments: 1) It's definitely not obvious that 4(a) has a Gaussian distribution. 2) Why would you expect 2a to be Gaussian?  The variance should be entirely due to the randomness in the selection of training data, which has nothing to do with Gaussians in the 0 noise case, as far as I can tell.
- In Figure 5, are the error bands the standard error of the mean, or the standard deviation across trials, or something else?
- I don't at all understand the setup of the 'Examining the neural network representation'.  I understand what the embedding of x is for the L2 model but what does that mean for the Implicit model, in which an x cannot be embedded without a y?