An implicit function learning approach for parametric modal regression
Authors: Yangchen Pan
Martha White
Amir-massoud Farahmand
Venue:  (2020)
Keywords: regression
modal regression
implicit function theorem
multivalue function

Abstract: 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 unimodal 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. 

Under review as a conference paper at ICLR 2020
AN IMPLICIT FUNCTION LEARNING APPROACH FOR PARAMETRIC MODAL REGRESSION
Anonymous authors Paper under double-blind review
ABSTRACT
For multi-valued functions such as when the conditional distribution on targets given the imputs is multi-modal-sttaaee standard regression approaches are not always desirable because they provide the conditional mean. Modal regression approaches aim 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 develop objective for learning 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.
INTRODUCTION
The goal in regression is to find the relationship between the input (observation) variable X E X and the output (response) that there exists an unknown underlying function g" X + y that maps the input space X' to the output space y. We only observe a noise-contaminated value of that function: sample (x, y) has = g" (x:) + for some noise 1j. If the goal is to minimize expected squared error, it is well known that IE] Y [x] is the optimal predictor (Bishop) 2006). is common to use Generalized Linear Models (Nelder & Wedderburn 1972). which attempt to estimate E| [x]] for different uni-modal distribution choices for p(y .). such as Gaussian (l2 regression) and Poisson (Poisson regression). For multi-modal distributions, however, predicting E[Y]z] may not be desirable, it may correspond to rarely observed y that simply fall between two modes. Further, this predictor does not provide any useful information about the multiple modes.
y variable, given samples of (X, Y). The underlying premise is
Modal regression is designed for this problem, and though not widely used in the general machine learning community, has been actively studied in statistics. Most of the methods are non-parametric. and assume a single mode jae Lee (1989); Lee & Kim [ Kemp & Silva (2012); Yu & Aristode- mou 2012); Yao & Li (014)) Lv et al. [eng Feng al. 2017). The basic idea adjust target values towards their closest empirical conditional modes, based on a kernel density estimator. These methods rely on the chosen kernel and may have issues scaling high-dimensional data due to issues in computing similarities in high-dimensional spaces. There is some recent work using quantile regression to estimate conditional modes (Ota al. 2018), and though promising for a parametric approach, is restricted to linear quantile regression.
parametric approach for modal regression would enable these estimators to benefit from the advances in learning functions with neural networks. The most straightforward way to do is to learn a mixture distribution, such as with conditional mixture models with parameters learning by a neural network (Powell) [oowell.198. 1987 Bishop 1994] Williimms [996] Husmeier 1997; Husmeier & Taylor 1998 Zen& & Senior 2014| Elleffsenn et al.) 2019). The conditional modes can typically be extracted from such models. Such strategy, however, might be trying to solve a harder problem than is strictly

---

Under review as a conference paper at ICLR 2020
needed. The actual goal is to simply identify the conditional modes, without accurately representing the full conditional distribution. Training procedures for the conditional distribution can be more complex. Methods like EM can be slow (Vlassis & Krose, 1999) and some approaches have opted to avoid this altogether by discretizing the target and learning a discrete distribution (Weigend & (Srivastava 1995; Feindi, (2004). Further, the mixture requires particular probabilistic choices to be made, including the number of components, which may not be correctly specified: they might be more or less than the true number of conditional modes.
In this paper, we propose a new parametric modal regression approach, by developing an objective to learn parameterized function fzz y) on both input feature and target/output. We use the Implicit Function Theorem (Munkres 1991), which states that if we know the onput- relation in the form of an implicit function, then a general multi-valued function, under certain gradient conditions, can locally be converted to a single-valued function. We learn function fa y) that approximates such local functions, by enforeing the gradient conditions. We empirically demonstrate that our method can effectively learning the conditional modes several synthetic problems, and that for those same problems, scales well when the input is made digh-diimnsional. We also show an interesting benefit that the joint representation learned over and y appears improve prediction performance even for uni-modal problem, for high frequency functions where the function values changes quickly between nearby æ. Finally, we show that our method provides small improvements that have asymmetric distributions over the targets. The proposed approach to multi-valued prediction flexible, allowing for a variable number of conditional modes to be discovered for each x, and we believe it is a promising direction for further improvements in parametric modal regression.
two regression datasets
2 PROBLEM SETTING
We consider a standard learning setting where we observe dataset of n samples, S { y/ }) Instead of the standard regression problem, however, we tackle the modal regression problem. The goal in modal regression is to find the set of conditional modes
M(x) y Öp(z. v) 0, p'p(z. v)
(1)
dy
Oy²
(x) is in general a multi-valued function. Consider the example in Figure [ There are two conditional modes for given æ. For example, for 0, the two conditional modes are yı -1.0 and 1/2 = 1.0.
The standard approaches to find these conditional modes involve learning p(y[z) to use non- parametric methods directly estimate these conditional modes. For example, for a conditional Gaussian Mixture Model, relatively effective approximation of these modes are the means of the conditional means. More generally. to get precise estimates, non- aanpaaammetric algorithms are used, like the mean-shift algorithm (Yizong Cheng 1995). We refer readers to Chen 2018); Chen et al. (2014) for a detailed nice review. These algorithms attempt to cluster points based on x and y, to find these conditional modes.
Looking at the plot in Figure I however, natural idea is to instead directly learn parameterized
function f(z, y) that captures the relationship between x and y. Unfortunately, it not obvious how to do so, nor how to use f(x,, y) to obtaín the conditional modes. In the next section, we develop an approach to learn a parameterized (xx,)) that can used to extract conditional modes, by using the implicit function theorem.
3 AN IMPLICIT FUNCTION LEARNING APPROACH
In this section, we develop an objective to facilitate learning parametric functions for modal regression (for value- prediction). The idea is to directly learn a parameterized function f(x y). where the set of minimum y for a function of f(z y) corresponds to the conditional modes. The approach allows for variable number of conditional modes for each r. Further, it allows us take advantage of general parametric function approximators, like neural networks, to identify these modal manifolds that capture smooth relationship between the conditional modes and
Consider learning an f(z y) such that f(z, y) = for all conditional modes and norzzero otherwise. For example, for the Circle problem, f(I,)) y) 2 y2 for all conditional modes y. Such
2

---

Under review as a conference paper at ICLR 2020
y
p(y|x=0)
0
0
X
(a) A dataset looks like circle
(b) Conditional distribution p(y/æ 0)
Figure 1: (a) shows dataset generated by uniformly sampling x E (-1, 1), conditioned which y is sampled from 0.5N (vI z², 00.... 0.11).. (b) shows the conditional distribution when æ = 0.
0.13).
strategy-finding fa, y) for all conditional modes-is flexible in that allows for different number of conditional modes for each æ. The difficulty with learning such an f. particularly under noisy data, is constraining to be zero for conditional modes yj and non-zero otherwise.
To obtain meaningful conditional modes y1..... Vm, for the y around each yj should be described by the same mapping (z). The existence of such 9j guaranteed by the Implicit Function Theorem (Munkres) 1991), under one condition on f.
Implicit Function Theorem: Let f Re IR. R be continuously differentiable function. Fix point (æ, y) E R4 x I&* such that f(x, y) = 0, for O E R If the Jacobian matrix J, where the element in the ith row and jth column is J|y] Byli invertible, then there exists open sets u, V containing (2,, y) such that there exists unique continuously differentiable function U + V satisfying g(x) = y and ((x, g())) O.
The theorem states that we know the relationship between independent variable æ and dependent variable y in the form of implicit function f(x,8/ y) O, then under certain conditions, we can guarantee the existence of some function defined locally to express given x. For example, a circle on two dimensional plane can be expressed {(x,, v)]=² + y²² 1}. but there no definite expression (single-valued function) for in terms of x. However, given a specific point on the circle (xo, 30) (other than Yo). there exists an explicit function defined locally around (In. Yo) express in terms of ær Notice that at y/ 0, the condition required by the implicit function theorem is not satisfied: 0(x + +--"-11
2y at y/ 0, and so is not invertible.
By
Obtaining such smooth local functions enables us to find these smooth modal manifolds. The conditional modes 91...., Im, satisfy f(x, 9 (x)) 0 and (f(x.g((*)) * 0. When training f. we can attempt satisfy both conditions to ensure existence of the Jj- The gradient condition ensures that for locally around (æ, 9 (a:),, we have (x, y) * 0. This encourages the other requirement that f(r, y) be non-zero for the that not conditional modes. Notice, though, that this condition only local and may not encourage a minimal set of conditional modes. We can also add a negative sampling component that encourages all y to default have f(x O. We find empirically, however, that this addition is not necessary. We therefore pursue this simpler objective, which avoids the complications of negative sampling.
0y
Now we derive the full objective, under stochastic targets. To do so, we make some assumptions the noise around the conditional modes. In particular, assume that the noise around each conditional mode is Gaussian. More precisely, define
(X,Y)
del (x(X)
3

---

Under review as a conference paper at ICLR 2020
for gj the conditional mode for Y. Our goal is to approximate c(x, y) with parameterized function fo(2,y) y) for parameters 8. We assume
<(X,))²)
<x, ~
= 0, ~
(2TT) k/2 exp
(3)
20
We also assume that dc(a.v)] for i + j: the conditional mode for dimension i is not influenced by グツ
changes to the in dimension j. Consequently, the implicit function condition becomes de(a.y(3))
A 0,
0bbl
to ensure invertibility. In fact, for this function, we have
Ōe(x, v) (g)(x) (z) - Oy
(4)
8y
Therefore, when learning f. we simply need to constrain (fl(2.9) -1 for all y. Putting this By
together, our goal is to minimize the negative log likelihood of fo. which approximates zero-mean Gaussian random variable, under this constraint which encourage with quadratic penalty term. This gives the following objective, where the goal is to find arg ming L(0):
@ S v.) dy
I(O) del
fo(
(5)
The same objective is used when doing prediction, but now optimizing over Given x²" we compute a y" in the set arg min, fo("" -))²+ modes, because the objective should be minimal for conditional modes. One option is to use gradient descent to find these y" In all our experiments, we opted for the simple strategy of searching over 200 evenly spaced values in the range of y.
(Qfu(e' + +1 These y" should correspond the conditional
0g
4 THE PROPERTIES OF IMPLICIT FUNCTION LEARNING
In this section, we conduct experiments to investigate the properties of our learning objective. First, by using the Circle datasets, we show its utility for dealing with multimodal distribution, particularly compared to modeling the entire distribution with mixture distributions. Second, inspired by the first experiment, we show that our algorithm achieves superior performance when the underlying true function is high frequency function. We empirically show that this is because our algorithm can better leverage the NÑ's representation power than the l2 regression. Hence, even with regular single-modal datasets, our algorithm could still be beneficial in some cases.
4.1 THE CIRCLE DATASETS
The purpose of this experiment is to study the properties of our algorithm for modal prediction problems. Mixture Density Networks (MDN) (Bishopl 1994) are used as baseline, compare to an approach that learns the distribution. For both algorithms, we use two hidden layer neural network (16 16 tanh units) and train by 228-mini-batch-sin stochastic gradient descent. We optimize both algorithms by sweeping learning rate from {0.1, 0.01, 0.001, 0.0001 For the purpose of evaluating prediction performance by testing error, we compute the root mean squared error (RMSE) for the predicted value and the true value, and the true value is defined as the one closer to the predicted value. For MDN to predict mode, given a point x, we search over 200 evenly spaced y values to maximize the learned log likelihood. We refer readers the Appendix B for study of the sensitivity of parameter in MDN.
We conduct two experiments datasets respectively. On the latter, we further conduct a challenging experiment by projecting the one dimensional feature value to 128 dimensional binary feature through tile coding" !I] We empirically show that: 1) our algorithm (Implicit) achieves higher sample efficiency than MDN; 2) on both datasets, our algorithm uses the same parameter settings, while MDN is quite sensitive to the number of mixture components and has a larger variance across different random seeds; 3) our algorithm can scale to high dimensional feature space and still maintain good performance.
single-circle (two modals) and double-circle (four modals most)
We refer readers to http://www..tt http://www.ncomppleo wwwpommpeeeeiias.eeeiieeeee net/:tless//ttm/// for more details about tile coding. It is frequently used feature generation method in reinforcement learning.
4

---

Under review as a conference paper at ICLR 2020
8.5
0.5
Implicit
Implicit MDN
Implicit
6.4
MDN
MDN
Noo suars
toxt
can uno
eersggad sorurs
avareedd
souns
0.1
6.0
5000
19000
15000
23000
0.0
0
5000
10d00
15000
20000
10000 20000 30000 40000 50000
(a) Single-cirele
(b) Double-cirele
(c) highdimensional doublc-circle
Figure 2: (a)(b) shows learning curves on single-circle and double-circle datasets respectively. In the single-circle problem, MDN uses 3 mixture components for learning; on the double-circle dataset and the high dimensional variant, uses 4 components to learn. All results are averaged over 30 random seeds and the shaded area indicates standard error.
The training set is acquired by uniformly sampling 40, 000 data points from circle: {(x,y))a² + y² 1). Note that, since we add zero mean Gaussian noise with standard deviation o= 0.1 to targets, we can interpret the conditional probability distribution p(y(æ)
Single-circle.
two-component mixture Gaussian 2())2=) 0.5. v1 x², o²) 0.5N -v1 xx2 o²) as shown in Figure On this dataset, though there are only two modals, MDN does badly (i.e. the testing error is almost outside of figure) with only two mixture components. hence
three components for it. The learning curves are shown in Figure 2 One can see that our algorithm learns much faster and is more stable in term of the standard deviation across different random seeds. Figure B (b) show the predictions outputted by our approach and MDN at the end of training (i.e. after 20000 updates).
Double-circle. On the double-circle dataset, the same number of training points are randomly sampled from two circles (i.e. {(z, y)k=² + y² 1}, (x, v)[=² + y² = -4}) and with the same Gaussian noise added to targets. This should be challenging dataset where p(y(z) can considered as piece-wise mixture of Gaussian: there are four components
E () and two components on x (-2, , -1) U 21, As a result, we set number components as 4 for our competitor MDN and keep all the parameter setting the same for algorithm. It should be noted that, the algorithm MDN learns much slower when moving from the two-components single-circle dataset to the cou-ncemso double-circle datasets. Additionally, components than the true distribution actually has to achieve superior performance, Similar to the single-cirele case, Figure3cc(d show the predictions the end of training, from which can see that the our algorithm does significantly better than the competitor.
find MDN requires more number of
High dimensional double circle. After using tile coding, we project the original one dimensional E [-2, 2] to binary features ((x) (0, 1) 11128 Then the projected features are used as input to the neural networks. One can see that although both algorithms converge faster, MDN converges to a worse solution than our algorithm does. To our best knowledge, this is the first time modal regression tested on such high dimensional dataset; most modal regression algorithms have been tested on low-dimensional problems.
Verify the learned error distribution. On both datasets, we investigate the validity of the as- sumption on the error distribution c(X, YY, by examining the empirical density of our learned error function fo(z, y) using all examples from training set. Figure 4 shows the eгтor distribution on the training set without noise (a) and with noise (b) added to the target variable. obvious that the error function indeed looks Gaussian distributed and the one trained without noise (4a)) shows an extremely small variance.
4.2 ROBUSTNESS TO HIGH FREQUENCY DATA
The above circle example can be thought of as an extreme case where the underlying true function has extremely high, even unbounded frequency (i.e. when the input changes little, there is
5

---

Under review as a conference paper at ICLR 2020
16
10
(a) Implicit, single circle (b) MDN, single circle (c) Implicit, double circle (d) MDN, double circle
Figure 3: (a)(b) shows the predictions of both our algorithm and MDN on the single circle dataset.
(c) show the predictions of the two algorithms on the double circle dataset.
120
4.0
3.5
100
3.0
80
2.5
60
2.0
1.5
40
1.0
20
0.0
-0.10
0.00
0.05
0.10
0.0
(a) No noise
(b) 0.1
Figure 4: (a) shows the empirical density of fo(x, y) for (Ir, y)s in the training set S when the algorithm is trained by not adding noise to target and adding noise (standard deviation o = 0.1)
respectively.
sharp change for the true target). As result, in this section, provide less extreme dataset and show that our algorithm has advantage dealing with such problems.
We generate synthetic dataset by uniformly sampling E [-2.5, 2.5] and using the below underlying procedure to compute the targets:
(i + x E [-2.5, + g x E 2.5]
sin(0.5mmz)
where & is zero-mean Gaussian noise with variance 2 This function has relatively high frequency when x E [-2.5, 0) and has relatively low frequency when x E [0, 2.5).
Significance of this artificial dataset. The dataset is designed to be difficult to learn, because several existing works (Smale et al.] 2004; [Smale & Zhou 2005 Jiang 2019) indicate that the bandwidth limit of the underlying true function strongly affects the sample efficiency of a learning algorithm. Intuition about why high frequency functions (large bandwidth limit) are difficult to learn can be gained from the Shannon sampling theorem, which characterizes the relation between minimum sampling rate to perfectly recover signal and the bandwidth limit of the signal (Zayed 1993): the sampling rate should exceed twice of the maximum frequency of the signal to guarantee perfect signal reconstruction.
Examining the learning behaviour. We use 16 x 16 hidden tanh units NN for our algorithm. For the l2 regression, we use the same size NN and perform extensive parameter sweep opti- mize its performance: activation function type swept over tanh and relu, learning rate swept over (0.1, 0.01, 0.001, 0.0001, 0.00001). For both algorithms, we used a mini-batch size of 128. Note the only difference between the two NNs is that Implicit has one more input unit, i.e. the target
In signal processing, sampling refers to the reduction of a continuous-time signal to discrete time signal. Sampling refers number of samples per second.

---

Under review as a conference paper at ICLR 2020
6.65
0.45
0.46
9.60
Implici
0.60
Implict
0.60
Implici
L2
0.40
0.40
marged
0.20
0.20
0.60
0
200560 socoss c60000 ascoco 1000000
0.00
0
200000
500009
650300 505000
1000000
0.60
200060
400009 660600 806000 1080060
(a) No noise
(b) a 0
(c) 0 = 0.2
Figure (a)(b)(c show performances of Implicit (red) and 12 regression L2 (black) objective as increase the Gaussian noise variance. We show the testing error measured by RMSE on entire testing set (solid line), on high frequency region (i.e. x E [-2.5, 0.0), dashed line) and on low frequency region (xr [0.0, 2.5]. dotted line). The results are averaged over 30 random seeds.
Figure 5:
Figure 5 (--)) shows the evaluation curve on testing set for the above two algorithms as the noise variance increases. The learning curve is plotted by evaluating the testing error every 10k number of iterations (i.c. mini-batch updates) averaged over 30 random seeds. We show the testing error on the entire testing set, on high frequency area (x S 0) and low frequency area (æ 2 0) respectively. We run into 1 million iterations to make sure each algorithm is sufficiently trained and both early and late learning behaviour can be examined.
Notice that, trained without observation noise (i.e. { 0), our implicit function learning approach achieves a much lower error (at the order of 10 than the la regression does (at the order of 10-1). As noise increases, the targets likely become less informative and hence algorithm's performance decreases to be closer to the l2 regression. Unsurprisingly, for both algorithms, the high frequency area is much more difficult to learn and is dominant source of the testing error. After sufficient training, our algorithm can finally reduce the error of both the high and low frequency regions to a similar level.
Examining the neural network representation. We further investigate the performance gain of our algorithm by examining the learned NN representation. We plot the predictions in figure a) and the corresponding learned NN representation through heatmap in figure In trained NN, we consider the output of the second hidden layer as the learned representation, and investigate its property by computing pairwise distances measured by l2 norm between 161 different evenly spaced points on the domain x E [-2.5, 2.5]. That is, point (r, x²) on the heatmap in figure b(b) denotes the corresponding distance measured by norm between the NN representations of the two points (hence the heatmap shows symmetric pattern w.r.t. the diagonal).
The representations provide some insight into why Implicit outperformed 12 regression. In figure 6 a), the l2 regression fails to learn one part of the space around the interval [-2.25, -1.1]. This corre- sponds to the black area in the left heatmap, implying that the l2 distance between NN representations among those points are almost zero. Additionally, one can see that the heatmap of our approach shows clearly high resolution on the high frequency area and low resolution on the low frequency area, which coincides with our intuition for a good representation: in the high frequency region, the target value would change a lot when changes a little, so we expect those points have finer representations than those in low frequency region. This experiment shows that given the same NN size, our algorithm better able to leverage the representation power of the NN.
RESULTS ON STANDARD REGRESSION DATASETS
In this section, we show that our approach can still be effective even for two standard real world datasets, with comparable performance to standard regression approaches. Our appendix B. includes learning curves and all details for reproducing the experiments. We compare to l2 and Poisson regression. For our algorithm, size NN but we consider hidden unit types as meta-parameter and we optimize them over tanh and relu. We report both root mean squared error (RMSE) and mean absolute error (MAF) on training
use 64 x 64 tanh units NN. For the la regression,
the same
7

---

Under review as a conference paper at ICLR 2020
0.5
0.0
-0.5
L2
Implicit
True
(a) Approximate and true functions
(b) Distance heatmap: L2(left) and Implici((right)
Figure 6: (a) Approximated functions and true function. (b) The distance matrix showed in heat map computed by hidden layer representation learned by L2 (left) and Implicit (right) method.
Table Prediction errors on bike sharing dataset. All numbers are multiplied by 102.
Algorithms
Train RMSE
Train MAE
Test RMSE
Test MAE
LinearReg
10094.40(+13.60)
7517.64 19.995
(0129.40(+:99.6))
7504.22(+ 44.20)
LinearPoisson
8798.26(+14.58)
5920.99(#13.66)
8864.90(=66.07)
5935.00(= 38.32)
NNPoisson
1620.46(147.71)
1071.39(129.55)
4150.03(177.76)
2616.49( 20.45)
L2
708. (08.2277)))79)
550.14(123.64))
3854.10(139.33
2560.83( 18.30)
Implicit
880.90(#30.53)
691.10(+23.99)
3683.76(157.12)
2426.52( 23.81)
and testing set respectively. All of our results are averaged over runs and each run, the data is randomly split into training and testing sets. Algorithms and datasets are as follows.
LinearReg. Linear regression, i.e. ordinary least square objective. The prediction is linear in term of input features. We use this algorithm reference.
LinearPoisson. The mean of Poisson is parameterized by linear function in term of input feature. NNPoisson. The mean of Poisson is parameterized by neural network (NN) (Fallah et al.) 2009).
Bike sharing dataset. We use the bike sharing dataset Fanaee-T & Gama (013) where the target Poisson distributed to show our algorithm's 's generality. Our prediction task to predict count of rental bikes in specific hour, given the information about weather, temperature, date, etc. We show the training result in table [ with standard errors. One can see that in the case of using linear function approximator, LinearPoisson has clear advantage over I LinearReg: however, in deep learning setting, we found NNPoisson achieves a worse performance than the regular la regression. Our algorithm, which does not make any assumptions the conditional distribution P(yla:), achieves slightly better performance.
Song year dataset. The song year dataset (Bertin-Mahieux al. 201 contains about half-million instances. The task is to predict a song's 's release year by using audio features of the song. The dataset has target distribution for which is not obvious which generalized linear model we should use. Hence we treat it as a regular regression dataset. One can see from Table that that our algorithm can slightly outperform regression, again potentially because we are not making distributional assumptions on p))).:)
Table 2: Prediction errors song year dataset. All numbers multiplied by 102.
Algorithms
Train RMSE
Train MAE
Test RMSE
Test MAE
LinearReg
681.56(+0.68) 95.340(=1.37)
157.56(+1..49)
681.66(+:1.55 1.52)
L2
798.61(+1.44))
563.51(10.65)
879.48(+1.744
606.57 1.995)
Implicit
822.53(+1.59)
582.11(+4.18)
869.62(+1.33)
600.73(42.15)
8

---

Under review as a conference paper at ICLR 2020
6 CONCLUSION AND DISCUSSION
The paper introduces a simple and powerful implicit function learning approach for modal regression. We show that it can handle datasets where the conditional distribution p(y(æ: is multimodal, and is particularly useful when the underlying true mapping has a large bandwidth limit. We also illustrate that our algorithm achieves competitive performance on large real world datasets with different underlying target distributions. We would like to conclude with the following future directions.
First, would be interesting to establish connections to KDE-based modal regression methods, which have a nice theoretical interpretation (Feng et al., 2017). The connection may yield finite sample analysis for our implicit function learning algorithm. Second, like many supervised learning algorithms, our algorithm may also overfit to noise. Popular regularization technique such random dropout (Srivastava al., (014) may be tested for very noisy data. Third, in online learning setting, the efficiency of doing prediction by arg min,, (1, y)² + (#f((z))) + 1)² becomes a concern. One possible solution is to borrow ideas from cross-entropy method as used in reinforcement learning (Lim) al. 2019). For example, we can use separate NN suggest
al. 2018; Simmons-Edle of initial values of for searching optimums by gradient methods. Last, it is worth investigating alternative constraints on the Jacobian instead of restricting the diagonal values 11
REFERENCES
Martín Abadi, Ashish Agarwal, Paul Barham, and et al. TensorFlow: Large-scale scale machine learning on heterogeneous systems, 2015. Software available from tensorflow.org.
Thierry Bertin-Mahieux, Daniel P.W. Ellis, Brian Whitman, and Paul Lamere. The million song dataset. In Proceedings of the 12th International Conference on Music Information Retrieval (ISMIR 2011), 2011.
Christopher M. Bishop. Mixture density networks. Technical Report NCRG/94/001, 1994.
Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
Yen-Chi Chen. Modal regression using kernel density estimation: A review. Wiley Interdisciplinary Reviews: Computational Statisties, 10(4):e143) 2018
Yen-Chi Chen, Christopher R. Genovese, Ryan Tibshirani, and Larry Wasserman. Nonparametric modal regression. arXiv e-prints, 2014.
Kai Olav Ellefsen, Charles Patrick Martin, and Jim Tørresen. How do mixture density mns predict the future? CORR, 2019.
Nader Fallah, Hong Gu, Kazem Mohammad, Seyyed Ali Seyyedsalehi, Keramat Nourijelyani, and Mohammad Reza Eshraghian. Nonlinear poisson regression using neural networks: simulation study. Neural Computing and Applications, 18(8):939, Jul 2009.
Hadi Fanaee-T and Joao Gama. Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, pp. 1-15, 2013.
Michael Feindt. A neural bayesian estimator for conditional probability densities. 2004.
Yunlong Feng, Jun Fan, and Johan A. K. Suykens. A statistical learning approach to modal regression. 2017.
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Yee Whye Teh and Mike Titterington (eds.), Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, volume 9 of Proceedings of Machine Learning Research, pp. 249-256 Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR.
Dirk Husmeier. Modelling conditional probability densities with neural networks. 01 1997.
Dirk Husmeier and John G. Taylor. Neural networks for predicting conditional probability densities: Improved training scheme combining and rvfl. Neural Networks, 11 11:89-116, 1998.

---

Under review as a conference paper at ICLR 2020
Myoung jae Lee. Mode regression. Journal of Econometrics, 42:337-349 1989.
Hui Jiang. A new perspective on machine learning: How to do perfect supervised learning. volume 0bss/900.000466 2019.
Gordon C.R. Kemp and J.M.C. Santos Silva. Regression towards the mode. Journal of Econometrics, 170(1):92 101, 2012.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. CoRR, abs/1412.6980, 2015.
M.-J. Lee and H. Kim. Semiparametric econometric estimators for truncated regression model: review with an extension. Statistica Neerlandica, 52(2):200-225, 1998.
Sungsu Lim, Ajin Joseph, Lei Lc, Yangchen Pan, and Martha White. Actor-expert: A framework for using action-value methods in continuous action spaces. CoRR, abs/18 10011/80009113 2018.
Zhike Lv, Huiming Zhu, and Keming Yu. Robust variable selection for nonlinear models with diverging number of parameters. Statistics Probability Letters, pp. 90-97, 2014.
James R. Munkres. Analysis On Manifolds. Boca Raton: CRC Press, 1991.
A. Nelder and R. W. M. Wedderburn. Generalized linear models. Journal of the Royal Statistical
Society, Series A (General), pp. 370-384, 1972.
Hirofumi Ota, Kengo Kato, and Satoshi Hara. Quantile regression approach to conditional mode estimation. arXiv e-prints, 2018
M. D. Powell. Algorithms for approximation. chapter Radial Basis Functions for Multivariable Interpolation: A Review, pp. 143-167. Clarendon Press, New York, NY, USA, 1987. ISBN 0-19-853612-7.
Riley Simmons Edler, Ben Eisner, Eric Mitchell, H. Sebastian Seung, and Daniel Lee. Q-learning for continuous actions with cross-entropy guided policies. CoRR, 1tt//1933.100000 2019.
Steve Smale and Ding-Xuan Zhou. Shannon sampling i:: Connections to learning theory. Applied and Computational Harmonic Analysis, 19(3):285-300 2005.
Steve Smale, René Thom, and Ding-Xuan Zhou. Shannon sampling and function reconstruction from point values. 2004.
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting. Journal of Machine Learning Research, 15:1929-1958. 2014.
N. Vlassis and B. Krose. Mixture conditional density estimation with the algorithm. In 1999 Ninth International Conference on Artificial Neural Networks ICANN 99, (Conf. Publ. No. 470), volume 2. pp. 821-825, 1999.
Andreas S. Weigend and Ashok N. Srivastava. Predicting conditional probability distributions: connectionist approach. International Journal of Neural Systems, 00:109-118, 1995.
Peter M. Williams. Using neural networks model conditional multivariate densities. Neural Computation, 8(4):843 8544, 1996.
Weixin Yao and Longhai Li. A new regression model: Modal linear regression. Scandinavian Journal of Statistics, 41 616666666,67, 2014.
Yizong Cheng. Mean shift, mode seeking, and clustering. IEEE Transactions on Pattern Analysis and Machine Intelligence, 17(8):790-777 799, 1995.
Keming Yu and Katerina Aristodemou. Bayesian Mode Regression. 2012.
A.I. Zayed. Advances Shannon's Sampling Theory. Taylor Francis, 1993.
H. Zen and A. Senior. Deep mixture density networks for acoustic modeling in statistical parametric speech synthesis. In 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp.
2014.
10

---

Under review as a conference paper at ICLR 2020
A APPENDIX
The appendix includes additional experimental results, and all experimental details for reproducible research.
B ADDITIONAL EXPERIMENTAL RESULTS
A classic inverse problem. One important type of applications of multi-value function prediction inverse problem. We now show additional results on a classical inverse function domain as used in (Bishop] 1994). The learning dataset is composed as following.
y 0.3 sin (2 + &, [0.1]
where is random variable representing noise with uniform distribution U(-0.1,0.1) We generate 80k training examples. In Figure 7 we plot the training dataset, and predictions by our implicit function learning algorithm with (arg min, fo y)² + @fa(a.v) + 1)²). We search over 200 evenly spaced ys in [0, 1] for 200 evenly spaced E [0, 1] get points in the form of (z, y)s.
By
0.6
0.0
0.0
0.2
0.4
0.6
0.8
1.0
(a) Training data
(b) Implicit
Figure 7: Figure (a) shows what the training data looks like. (b) shows the predictions our implicit leamning approach.
Additional result of mixture density network. We now provide additional result of mixture density network MND on the three multivalue prediction dataset:circle, double circle, and high dimensional double circle. The learning curves are shown in Figure g For each component on each training dataset, we choose stepsize from {0.1, 0.01, 0.001, 0.0001
MDN-2
MDN-3
mean
MDN-4
monn fmr
moat
MDN-6
0.1
0.1
0.1
5000
10006
15000
20000
5000
10000
15000
20000
10000 20000 30000 43000
(a) single circle
(b) double circle
(c) high dimensional double circle
Figure We show the learning curves of mixture density network with different number of components. In (a) and (b), the learning curves of component 2 out of figure. For each number of component, we choose the learning curve by optimizing other parameter (i.e. stepsize). The results averaged over 30 random seeds.
11

---

for Under review as a conference paper at ICLR 2020
B.1 REPRODUCE EXPERIMENTS IN THE PAPER
In this section, we provide additional information about datasets we used and experimental details for reproducing results in this paper. Our implementation based on Python 3.3.6. Our deep learning implementation is based on Tensorflow 1.11.0 et al. (2015). All of our algorithms trained by Adam optimizer (Kingma & Ba. 2015 with mini-batch size 128 and all neural networks initialized by Xavier (Glorot & Bengio 2010 For implicit function learning algorithm,
use tanh units for all nodes in neural network We search over 200 evenly spaced values for prediction
except song year dataset where use 100. Best parameter settings used to reproduce experiment
on each dataset are showed in figure 9 the end of learning.
The best parameters are chosen according to the testing error
Algorithms datasets
sharing
Song year
frequency test, eq(6)
circle/double circle
problem
L2 regression
learning rate 0,0001. 64-by-64
learning rate
learning rate=0.01 16-by-16 tanh units
0.0001. 64-by-64
tanh units
tanh units
NN
Implicit
leamning 0.0001. 64-by-64
learning
learning
learning 16-by-16 tanh units NN
learning
0.0001, 64-by-64 tanh units
0.001, 16-by-16
0.01, 16-by-16 6 tanh units NN
tanh units
tanh units NN
LinearReg
learning 0.0001
learning rate 0.0001
PolssonReg
learning
NNPoisson
learning rate
0.0001, 64-by-64 tanh units, linear output unit
Figure 9: Best parameter setting for reproducing experiments.
B.1.1 CIRCLE AND DOUBLE CIRCLE EXPERIMENT
Circle dataset is generated by uniformly sampling x E [-1,1] first and then = v1 x² or y = - --1 2 with equal probability. Double circle dataset generated by uniformly sampling an angle or r = 2.0 with equal probability. High dimensional double dataset is generated by mapping the http://www.incommlleeeeeee.
[0, 2m] then polar expression to compute
a, y sin where 1.0
original x to (0, 128 dimensional space. We refer tile coding software. The setting of tile coding we used to generate feature is: memory size 128, 8 tiles and tilings. We keep the neural network size the same we used in low dimensional case, i.c. 16 x 16 tanh units.
For mixture density network (MDN), we use tanh hidden layers and three mixture components on single circle examples, four mixture components on double-circk example. We sweep over learning rate from (0.1. 0.01, 0.001, 0.0001) and the best it chooses 0.001 on the single circle dataset and 0.01 the double circle dataset. The maximization is done by using MLE, the method described in the original paper Bishop (1994).
1.1.2 HIGH FREQUENCY DATA EXPERIMENT
The dataset is generated by uniformly sampling x E [-2.5, 2.5) and then compute targets according to the equation 6
sin (8TZ)
x E [-2.5,0 x E [0, 2.5]
y =
sin (0.5mzz
'Rigorously, to satisfy the assumption that (x(,, Y) is Gaussian distributed, linear output unit should be used. However, we observe large error rarely happens. In fact, in our case, it is easy to see that assuming the distribution be truncated Gaussian (then using tanh is justified) would yield the same optimization objective.
12

---

Under review as a conference paper at ICLR 2020
Table 3: Prediction errors song year dataset with author suggested train-test split. All numbers are multiplied by 10a. The randomness comes from neural network initialization and stochastic mini-batch update.
Algorithms
Train RMSE
Train MAE
Test RMSE
Test MAE
LinearReg
956,43(=0005)
952.22(#00.9) 18).66(=0.31)
681.22(+0.39)
L2
836.77(=1047))
584.00(=1.22) 85.548(=0.51)
888.58(10.65) 886.03(#0.63) 610.52(+1.10)
612.79(10.67)
Implicit
593.04(#1.11)
We sweep over f0.1, 0.01, 0.001, 0.0001, 0.00001 } to optimize stepsize for both the 12 regression and our algorithm, while we additionally sweep over hidden unit and output unit type for the la regression from (tanh, relu). The best parameter is chosen according the testing error at the end of training, and the testing error is averaged over 30 runs and at each run, the data is randomly split into training and testing sets. Since the best learning rate chosen chosen by the l2 regression is 0.01 while our algorithm chooses 0.001, in figure I0 we also plot the learning curve with learning rate 0.001 to make sure that the performance difference
due to a slower learning rate of our algorithm.
0.65
0.60
Implicit
Implicit
L2
0.50
0.50
0.40
Bog. saave
0.30
0.20
0.10
0.10
1000000
400000
(a) la with learning rate 0.01
(b) La with learning rate 0.001
Figure 10: In (a) we repeat the figure shown in previous
B.1.3 EXPERIMENTS ON REAL WORLD DATASETS
The bike sharing dataset (Fanaee- T & Gama, 2013)
tdatase+d and song year dataset Bertin- Mahieux et al. (2011) (https://archive. www.ivviis..iiieuuu. cllaastt/yearprediiiionnn very different target distributions as shown in figure T2 On the two dataset, we use 64 64 tanh hidden units and sweep over learning rate from (0.01, 0.001, 0.0001 For the l2 regression, we found that using tanh unit works better than relu.
information are presented in figure [1] Note that the two datasets have
Notice that the contributor of song year dataset suggests using the last 51630 testing set (Berting Mahieux et al. (20TT) hence we also report this additional result in table 3 The performance actually quite similar to that by random split.
13

---

Under review as a conference paper at ICLR 2020
Datasel and preprocessing information
Number of instances
Train size
Text size
Inpul feature dimension after
Input feature preprocess
Target preprocess
preprocesingg
sharing
17379
remove
[0, except for poisson regression algorithms; scaled back when
date, index.
weather situation
weekday registered, casual; use osc-hot encoding for categorical variables
compute test error
Song Year
515345
412276 463715
103069 or 51630
90 standardize to unit variance: statistics acquired by using training set
Scale to [0, I]; scaled back when compute test erroe
Figure 11: Data preprocessing information.
0.175
0.012
0.150
0.010
0.125
0.008
0.004
0.050
0.002
0.025
0.000
0.000
200
400
600
800
1000
1920
1940
1960
1980
2000
(a) Bike sharing dataset target distribution
(b) Song year datasel target distribution
Figure 12: Bike sharing targets show clear Poisson distribution while song year dataset's target distribution is not intuitive.
14


Reviews:
[{"review_id": "Sye133lTKS", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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.\n\nThe 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. \n\nI 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.  \n\nMy 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.  \n\nI 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.\n\nSpecific questions:\n- In section 3, second sentence, I think you want to write either the minimum of f^2 or a zero of f?\n- Page 3, not sure what is meant by \"should be described by the same mapping g_j(x).\"\n- In the implicit function theorem, it ends with a capital F, but I think you want f.\n- 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?\n- In Equation 5, you write a full l_2 norm for what I believe is just a scalar -- how about parenthesis instead?\n- 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?\n- 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.  \n- 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.\n- In Figure 5, are the error bands the standard error of the mean, or the standard deviation across trials, or something else?\n- 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? \n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "N/A"}, "extra_text": {}}, {"review_id": "HklWwuchKB", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes a parametric modal regression algorithm for multi-modal distributions. In such settings, learning the conditional mode is more desirable than the conditional mean, the latter being the focus of standard regression approaches. The objective in model regression is to find the conditional mode. \n\nExisting non-parametric approaches to learning the conditional mode are difficult to scale. On the other hand, existing parametric approaches for modal regression aim to learn the full conditional distribution, which is a much harder problem. The present work proposes a parametric approach to estimate the conditional mode using the Implicit Function Theorem. The key idea is to learn a (parametric and implicit) function f(x,y) whose minima(s) 'y' corresponds to the conditional modes. However, the paper does not explain clearly how the above idea is put to practice as an algorithm in Section 3. It will be great if Section 3 can be re-written with more detailed explanations in paragraphs following the \"Implicit Function Theorem\" definition and more examples. \n\nThe experiments mainly are performed artificial datasets to assess the algorithm when the underlying multi-model distributions are known. Empirical results are also discussed on a couple of real-world datasets where the proposed algorithm obtains the best results. ", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_thoroughness_in_paper_reading": "I made a quick assessment of this paper.", "review_assessment:_checking_correctness_of_experiments": "I did not assess the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "HJlJbLHVtB", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 1, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper proposes an implicit function approach to learning the modes of multimodal regression. The basic idea is interesting, and is clearly related to density estimation, which the paper should have discussed. On the other hand the tackled problem is a weaker variant of the learning multimodal output distributions in regression. What is the benefit of tackling this simpler problem, when in reality we are interested in the output distributions? In Bayesian literature (BNNs, deep GPs, Bayesian regression) multimodal outputs are routinely produced, but the paper ignores these methods. The paper should compare to these, and demonstrate the benefit of just learning the modes instead of learning the full output distributions (perhaps the mode learning gives more accurate mode assignment?).\n\nThe paper should be more rigorous on what is the problem setting. Right now its unclear if g(.) is multimodal, if \\eta is multimodal, or if both are. This is an important distinction since it deals with whether the underlying system is modal or if the noise modal, and connects with aleatoric/epistemic characterisations. The earlier works could have been presented in more clear manner, currently the paper gives an impression that learning joint function over (x,y) is novel which it surely isn't. It was a bit hard to follow the earlier works, they should be more clearly categorised.\n\nThe artificial datasets raise more questions than they answer. Both of them are simple univariate cases, where one would expect even simple approaches work very well. But there are suddenly 40000 datapoints for the first case (why not eg 50?), optimisation seems to take hundreds of thousands of iterations, and the results are still quite bad all around (fig3). I suspect that there are coding issues, and overall the univariate examples are not very informative in general. The \"high-dimensional\" case is not really high-dimensional, and a real one should have been proposed instead. There are also no real competing methods in any experiments, except for very simple regression methods or baselines from the 90's (MDN). With no comparisons to competing methods the experiments are close to worthless. The two real dataset experiments are not enough, and even there the improvement seems to be modest with only trivial baselines included. I also don't understand why the bike classification problem is chosen for regression.\n\nThe paper proposes to combine the implicit function idea with output mode estimation, however the problem definition is vague, competing Bayesian methods and density estimation methods are ignored, the experiments are insufficient with little state-of-the-art comparisons, few datasets and clear problems in the learning. The results show only small improvements, which are not explicated sufficiently. The specific problem tackled is of minor importance for the wider ICLR community.", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_thoroughness_in_paper_reading": "I made a quick assessment of this paper."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes an implicit function approach to learning the modes of multimodal regression. The basic idea is interesting, and is clearly related to density estimation, which the paper does not discuss. 

Based on the reviews and the fact that the authors did not submit a helpful rebuttal, I recommend rejection.
Decision: Reject