You're really asking the wrong question here...
Let's back up a bit. You're attempting to estimate some parameters here, either by a maximum likelihood method or more likely by $\chi^2$ minimization, right? Normally, after you've completed the $\chi^2$ minimization, you take the Hessian of $\chi^2$ at the optimal parameters, and then use its inverse as a covariance matrix for the fitted parameters, right? However, you've found that your Hessian matrix is essentially singular, right?
So what does this mean?
First, the fact that the Hessian is essentially singular tells you that the $\chi^2$ minimization problem is ill-conditioned in the sense that many different sets of parameters fit the data equally well. This presents problems for whatever algorithm you used to perform the $\chi^2$ minimization, and it's likely in practice that you haven't got a very good minimum.
More importantly, the fact that the $\chi^2$ minimization problem is ill-conditioned tells you that the experiment you conducted really doesn't pin down the values of the parameters that you're trying to estimate.
It's important to understand that the problem you have here is not one of numerical computing, but that it is in fact a significant statistical issue. What you really need to do here is redesign your experiment so that it will provide you with useful information about the parameters that you're trying to estimate.
Trying to kludge together some sort of approximate covariance matrix by using the truncated SVD is simply the wrong way to approach this problem. What will go wrong? By truncating the SVD, you'll end up with confidence intervals for the fitted parameters that are overly tight.
Best Answer
Cholesky sounds like a good option for the following reason: once you've done the Cholesky factorization to get $C=LL^T$, where $L$ is triangular, $x^TC^{-1}x = ||L^{-1}x||^2$, and $L^{-1}x$ is easy to compute because it's a triangular system. The downsides to this are that even if $C$ is sparse, $L$ is probably dense, and also that you do the same amount of work for all $C$ and $x$ while other methods may allow you to exploit some special structure and get good approximations to the solution with less work. For those reasons, you might also consider Krylov subspace methods for computing $C^{-1}x$, like conjugate gradients (since $C$ is symmetric and positive definite), especially if $C$ is sparse. $n=250$ isn't terribly large, but still large enough that Krylov subspace methods could pay off if $C$ is sufficiently sparse. (There might actually be special methods for computing $x^TC^{-1}x$ itself as opposed to $C^{-1}x$, but I don't know of any.)
Edit: Since you care about stability, let me address that: Cholesky is pretty stable, as you note. Conjugate gradients is notoriously *un*stable, but it tends to work anyway, apparently.