Of course, my response should be - picky, picky, picky!
This is a consequence of a log transform I guess. You could use the 'springs' method as a regularizer to avoid that, but it might be too extreme.
So an alternative is to use a piecewise transformation. Thus, choose some breakpoint, above which the transformation will be linear, and below which, it will behave as a log transform. I'll call that value z0. So we can write the transform function W(z) as
W = @(z,z0) log(z).*(z<=z0) + ((z - z0)/z0 + log(z0)).*(z>z0);
You can plot it, picking perhaps z0=3 as a breakpoint for an example.
Of course, the transformation is invertible. This should work:
Winv = @(w,z0) exp(w).*(w<=log(z0)) + z0*(w - log(z0) + 1).*(w>log(z0));
Again, we can plot it, here for z0 = 3 again, as:
Something like that should give you the best of all worlds, and all you need to do is choose some value for z0 that makes sense to you. I might start with the mean of z as a good choice for z0.
Hi John. Thanks for the tip regarding using a log transform. It seems to be a sensible approach and I tried it. The fitted surface is now all positive, as expected. However, it also overshoots: the highest peak (Z) in the fitted surface is now 3 times the largest data point. Previously the fit closely followed the data.
I should add that the data did have 0's in it and I didn't want to completely remove them because they represent boundary conditions for the somewhat sparse data I am fitting to; for example, the spectrum vanishes at a certain energy (e.g. E=0) or at some angles. Instead, as you suggested, I replaced the 0's with a small number; I tried 1E-10, 1E-5, 1E-1, and 1. The overshooting appeared in all cases and is probably unrelated to the presence of 0's or small numbers in the data.
I tried to attach pictures of the fitted surface before and after the log transform, but there doesn't seem to be a way to do that on this forum.
Neil - the answer is definitely yes, and no. Or maybe it is the other way around.
No, it would be difficult to make gridfit work as you would like to see via internal changes, since that would be a bound constraint on a fairly massively large scale linear algebra problem. So while I could in theory use a tool like lsqnonneg or lsqlin instead of backslash for the solve, it would take forever to terminate on any problem of reasonable size.
Regardless, there is a solution that will work for you, and is quite simple. Typically when there is a bound constraint at zero because of physical issues like this it is because the problem really should be transformed to make it more linear. The logical transformation is a log transform. Thus instead of fitting Z as a function of X and Y, fit the surface in the form of log(Z) as a function of X and Y. Clearly there will no longer be any negativity issues, because to recover the surface you really want, you will simply exponentiate.
If you actually had any data points that were an exact zero, it would be best to drop them completely. Alternatively you could replace them with some small number, but that would cause the problem to be less smooth and might introduce bugs in the surface. Of course, the non-smooth region will be in areas where the log will be very negative, so after exponentiation it will be essentially zero anyway. So either way will work on those exact zeros.
I hope the log transform idea solves your problem. It is an approach I have used often with good success.
Thank you for making this code available. Is there a way to tweak it to add the additional constraint that the fitted surface should also always be positive? The X,Y,Z data I'm fitting represent spectrum (Z) vs. energy (Y) and angle (X). Therefore, negative Z is unphysical.