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. The goal in regression is to find the relationship between the input (observation) variable X ∈ X and the output (response) Y ∈ Y variable, given samples of (X, Y ). The underlying premise is 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 y = g * (x) + η for some noise η. If the goal is to minimize expected squared error, it is well known that E[Y |x] is the optimal predictor (Bishop, 2006) . It is common to use Generalized Linear Models (Nelder & Wedderburn, 1972) , which attempt to estimate E[Y |x] for different uni-modal distribution choices for p(y|x), such as Gaussian (l 2 regression) and Poisson (Poisson regression). For multi-modal distributions, however, predicting E[Y |x] may not be desirable, as 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. 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 (1998) ; Kemp & Silva (2012) ; Yu & Aristodemou (2012) ; Yao & Li (2014) ; Lv et al. (2014) ; Feng et al. (2017) . The basic idea is to 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 to 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 et al., 2018) , and though promising for a parametric approach, is restricted to linear quantile regression. A 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 so is to learn a mixture distribution, such as with conditional mixture models with parameters learning by a neural network (Powell, 1987; Bishop, 1994; Williams, 1996; Husmeier, 1997; Husmeier & Taylor, 1998; Zen & Senior, 2014; Ellefsen et al., 2019) . The conditional modes can typically be extracted from such models. Such a strategy, however, might be trying to solve a harder problem than is strictly 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; Feindt, 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 a parameterized function f (x, y) on both input feature and target/output. We use the Implicit Function Theorem (Munkres, 1991) , which states that if we know the input-output 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 a function f (x, y) that approximates such local functions, by enforcing the gradient conditions. We empirically demonstrate that our method can effectively learning the conditional modes on several synthetic problems, and that for those same problems, scales well when the input is made high-dimensional. We also show an interesting benefit that the joint representation learned over x and y appears to improve prediction performance even for uni-modal problem, for high frequency functions where the function values changes quickly between nearby x. Finally, we show that our method provides small improvements on two regression datasets that have asymmetric distributions over the targets. The proposed approach to multi-valued prediction is 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. 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|x) 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, it 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 as random dropout (Srivastava et al., 2014) may be tested for very noisy data. Third, in online learning setting, the efficiency of doing prediction by arg min y f θ (x, y) 2 + ( ∂f θ (x,y) ∂y + 1) 2 becomes a concern. One possible solution is to borrow ideas from cross-entropy method as used in reinforcement learning (Lim et al., 2018; Simmons-Edler et al., 2019) . For example, we can use a separate NN to suggest a set of initial values of y for searching optimums by gradient methods. Last, it is worth investigating alternative constraints on the Jacobian instead of restricting the diagonal values to −1.
We introduce a simple and novel modal regression algorithm which is easy to scale to large problems.  The paper proposes an implicit function approach to learning the modes of multimodal regression. The present work proposes a parametric approach to estimate the conditional mode using the Implicit Function Theorem for multi-modal distributions. 