I have a real function $f:\mathbb{{R}}^{d}\mapsto\mathbb{R}$, where $d>1$. The question is how to compute the level set $A=\left\{ \theta:f\left(\theta\right)\geq a\right\} $. I am a statistician knowing very few things about computer science. I would appreciate very much if anyone could suggest some reference, methods and/or alogirthms.

In my problem, $f\left(\theta\right)=\sum_{i=1}^{N}\log g\left(\mathbf{z}_{i};\theta\right)$, where $g\left(\mathbf{z}_{i};\theta\right)$ is a probability density function of $\mathbf{z}_{i}$ parameterized by a vector of parameters $\theta$. $\mathbf{z}_1,\ldots, \mathbf{z}_N$ are observed samples. For example, $g\left(\mathbf{z}_{i};\theta\right)$ is the density function of a normal distribution with variance $1$ $$ g\left(\mathbf{z}_{i};\theta\right)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\left(y_{i}-\mathbf{x}_{i}'\theta\right)^{2}}{2}\right). $$ Here $\mathbf{z}_{i}=\left(y_{i},\mathbf{x}_{i}\right)'$.

  • $\begingroup$ How is $f$ represented, i.e., what do you know about $f$? Clearly, there are uncountably many $f$s and only countably many methods/algorithms, so we need to know something about $f$. $\endgroup$
    – Pål GD
    Jul 20, 2013 at 19:33
  • $\begingroup$ Thanks for reply. I have updated my question with more details about the function $f$. Hopefully, the question is better framed. $\endgroup$
    – semibruin
    Jul 20, 2013 at 19:50
  • $\begingroup$ Real numbers can not be handled in general; what are your requirements on the result (quality)? $\endgroup$
    – Raphael
    Jul 20, 2013 at 20:51
  • 3
    $\begingroup$ @Raphael You are being very pedantic. Real numbers are manipulated all the time in computers, usually using a floating point representation. Of course, computations are approximate. It could be that semibruin is interested in solving a practical problem rather than some theoretical question. $\endgroup$ Jul 20, 2013 at 21:06
  • $\begingroup$ @YuvalFilmus: See, that's why I'm asking. I fail to see how it's "pedantic" to inquire what the desired result is, if only to rule out misconceptions. (And no, IEEE floats are not reals. And that really hurts in any kind of numerical computation.) $\endgroup$
    – Raphael
    Jul 21, 2013 at 17:20

1 Answer 1


We can compute $$ f(\theta) = -N\log (2\pi)-\frac{1}{2}\sum_{i=1}^N (\langle x_i,\theta \rangle -y_i)^2. $$ Expanding this out, we get that $f(\theta)$ is some quadratic form: $$ f(\theta) = \theta' A \theta + v'\theta + C, $$ where $A$ is symmetric. The next step is to get rid of the linear term. Let $\theta = \psi + \epsilon$. Then $$ \theta' A \theta + \langle v,\theta \rangle + C = \psi' A \psi + (2\epsilon' A + v') \psi + (\epsilon' A \epsilon + v'\epsilon + C). $$ So choosing $\epsilon' = -v' A^{-1}/2$ (assuming $A$ is invertible), we get $$ f(\psi + \epsilon) = \psi' A \psi + C'. $$ The theory of quadratic forms tells us that if $A$ is invertible then we can write $A = P'DP$, where $D$ is a diagonal matrix with $\pm 1$ on the diagonal and $P$ is invertible. If $\psi = P \chi$ then $$ f(P^{-1} \chi + \epsilon) = \chi' D \chi + C' = \sum_{i=1}^k \delta_i \chi_i^2 + C', $$ where $\delta_i \in \{\pm 1\}$. Suppose for example that the first $k_+$ are $+1$ and the rest $k_- = k - k_+$ are $-1$, and write $\chi = [\chi_+\;\chi_-]$. Then $$ f(P^{-1} \chi + \epsilon) = \|\chi_+\|^2 - \|\chi_-\|^2 + C'. $$ If $k_+ = 0$ or $k_- = 0$ then the level sets are spheres or cospheres, up to an affine transformation of course. Otherwise it gets more complicated.

In your case, however, it is clear from the function $f$ itself then the quadratic form $A$ is negative semidefinite, probably in fact negative definite. In other words, $k_+ = 0$. So the sets you're after are affine transformations of cospheres. In other words, they are "outsides" of ellipsoids.

The only non-trivial step is finding the matrix $P$. Since your matrix is negative definite, you can use Cholesky decomposition.

  • $\begingroup$ Thanks for your answer, Yuval. The normal distribution is an elliptical distribution. I guess your approach after slight modification shall work for all elliptical distributions. I am indeed working on more complicated scenario. Would it make some sense by considering the Taylor expansion $f(\theta)\approx\theta'f_{\theta=0}^{(2)}\theta+f_{\theta=0}^{(1)}\theta$. The hessian matrix here is usually negative definite in my case. $\endgroup$
    – semibruin
    Jul 21, 2013 at 6:17
  • $\begingroup$ The quality of the approximation depends on the problem, and whether it is good enough or not depends on the application. You'll just have to try it out. $\endgroup$ Jul 21, 2013 at 7:56
  • $\begingroup$ @semibruin Before implementing any formula/algorithm, make sure you understand how loss of precision propagates and what the precision of your final result is. That's not trivial, even if Yuval seems to think such considerations to be pedantic. $\endgroup$
    – Raphael
    Jul 21, 2013 at 17:22
  • $\begingroup$ @Raphael Thanks for your remark. I need to discuss the approximation error to make the method well established. Do you know any introductory or overview paper about determining precision? $\endgroup$
    – semibruin
    Jul 21, 2013 at 18:13
  • 1
    $\begingroup$ @semibruin Numerical stability is a delicate subject, and you're better off finding a library routine that solves the problem for you. Inverting matrices, for one, is looking for trouble. $\endgroup$ Jul 21, 2013 at 18:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.