The focus of the work is modeling non-negative functions (or more generally functions with output in convex cones). The authors propose a non-parametric approach to tackle this task. Particularly, they approximate non-negative functions by quadratic forms in reproducing kernel Hilbert spaces parameterized by positive semi-definite operators. They show that the proposed approach have various favorable properties including convexity, universal approximation, representer theorem (hence it is computationally tractable) and it gives rise to complexity bounds alike to standard kernel approaches. The efficiency of the approach is illustrated in density estimation, heteroscedastic Gaussian process estimation and multiple quantile regression. This a nice submission: the paper is clearly written, it delivers both interesting theoretical insights and have practical relevance.