Outstanding issues in greenspline

[PaulWessel]

Instead of making this an issue, I wish to tabulate mistakes and/or omissions in greenspline that we will need to address or at least document the limitations. Some of these have been fixed already but there are more.

Our example gspline_3.sh shows gridding of data, then calculation of the gradient of the solution in the az = 135 direction via -Q135. While the grey-shaded image is qualitatively correct (positive and negative in the right places and right direction of trend) the values are incorrect. This is because greenspline applies several adjustments and corrections to stabilise the inversion:

1. For many cases, including Cartesian 2-D, we fit a straight plane to the data since that shape is not reflected by the Green's functions. So we fit it, remove the trend, and operate on the residuals.
2. The residuals have a range min to max, and we normalise these residuals further by dividing by this range.

On output, we then evaluate the Green's function solution based off these residuals, then undo those corrections in the reverse order. This works fine. However, until very recently, we did not make the same corrections on input and output when we have gradient constraints or gradient output requests. Except for the 1-D spline case, which now does the right things, we still have work to do:

1. If there is a mixed data and gradient input data set, then we also need to determine the gradient of that plane (call it n_hat g, with n_hat the gradient direction unit vector and g the magnitude). Thus, at each input gradient point (with direction n_hat_j and magnitude g_j) we need to remove n_hat dot n_hat_j times g. In 1-D there is no direction so we just remove the slope of the LS line but in 2-D the directions matter, obviously.
2. If -Q is used then we internally evaluate the gradient of the solution but in the -Q direction. Still, we need to make adjustments for the sloping plane. I think if q_hat is the unit vector in the -Q direction we need to add back in q_hat dot n_hat times g to allow for the missing (plane) contribution in that forced direction. This is certainly not done.

So, master currently gives a different result for the spline_3.sh than the truth since the wrong undo normalisation takes place. The new result is closer to the correct answer but because the effect of the LS plane has not been considered it is not completed. Thus, I guess I could let the test fail (if we merge my current fixes) or I wait until I can do the q_hat stuff.

Another issue is what started this reexamination of greenspline in the first place: A user had both data and slope constraints at the same point - whereas we alway assumed the slopes were at different locations than data observations. Thus, in the current implementation, if -A is used these add m new rows in the Ax = b equation where the weighted sum of the gradients of the Green's functions times unknown coefficients match the observed slope. So each such point adds one more unknown. Yet, if slope and value are known at the same point then there is only one Green's function associated with that point - we just happen to have 2 constraints. This means we will have an overdetermined system. If there were 100 data points and 100 slope constraints at the same 100 points then we end up with 200 equations in 100 unknowns and need a least-squares solution (like we do if there are weights involved). This scenario is not implemented. I am not sure how to easily do this in terms of how we organize input (yet another input file for sloes that coincide in location with data?) or do we scan the -A file to determine if some (or all) points overlap exactly with locations in the primary input? Certainly, that would make it simpler for the user. However, a user may reasonably expect to give a single input file with both data and slopes on each row. So this needs more discussion I think. Also I do not have a need for this so there is a limit to how much time to spend on other peoples problems.

In terms of documenting limitations, here is a list:

• The -A option expects slope constraint not to coincide with data constraints in location.
• The -A currently is only known to work (i.e., tested) for 1-D cubic spline (-Sc)

[leouieda]

...regarding the remove-compute-restore of plane in gradients

Yes, agree 100% with the proposal.

if -A is used these add m new rows in the Ax = b equation where the weighted sum of the gradients of the Green's functions times unknown coefficients match the observed slope. So each such point adds one more unknown. Yet, if slope and value are known at the same point then there is only one Green's function associated with that point - we just happen to have 2 constraints.

A way to have both cases work:

1. Treat the gradient as new data, just like is currently being done.
2. But, instead of setting 1 node per data point, place nodes on the blockmean the data coordinates (values + slopes) at the desired grid resolution. Do this just for the node coordinates, not for the data values themselves. The data values stay the same.
3. If there are more data than nodes, solve with least-squares.

This way, if the data and slope points overlap (exactly or close-enough given spacing), you get the overdetermined system that is solved by least squares and save quite a lot of computation. This works since the least-squares will fit slope and value approximately and not exactly. If they don't overlap at all, you get the current behavior where you still get 1 node per data and slope.

In my experience, you really only need 1 node per final grid point. So this approach would even make the code faster for people who forget to do a blockmean before running greenspline.