I am fitting a 2D depth surface (z_mod) to some 170 data points (z_obs) using greenspline. To account for noise in the data, I’m using singular value decomposition to model the surface and am currently choosing a eigenvalue cut-off subjectively, using:
I’d now like to arrive at a more objective cut-off value by plotting the famous “L-curve” (data variance vs. model variance), frequently used for determining the damping parameter in inverse problems:
data variance
^
| | function of svd cut-off value
| |
| |
| |
| |
| \
| \* Optimal damping
| \_________________________
|
+--------------------------------------- >
model variance
I can compute the data variance as sum((z_mod - z_obs)**2). However, I’m having problem plotting the model variance. Essentially, I’m unable to access the spline coefficients alpha in the greenspline source.
I’ve tried various combinations of -E-C+f, but was unable to access the model parameters. Or am I overseeing something? Can I compute it from the eigenvalues (supplied by the -E option)?
Best fix: Is it me, have I missed something?
Dirty fix: Can I compile greenspline easily, so I can send alpha to stdout?
Elegant fix: Would it be possible to implement e.g. an option -C+afile: “For each cutoff eigenvalue, write out the alpha coefficients to file (we automatically insert “###” before the extension, see +c and +i).”
Dear community
I could compile the greenspline source without major issues and now have my private version in which I write out the model norm as an INFORMATION.
for (e = 0; e < (uint64_t)n_use; e++) { /* Only loop over the first n_use eigenvalues (if restricted) */
GMT_Report (API, GMT_MSG_INFORMATION, "Evaluate spline for eigenvalue # %d\n", (int)e);
gmt_M_memcpy (s, ssave, nm, double); /* Restore original values before call */
(void)gmt_solve_svd (GMT, A, (unsigned int)nm, (unsigned int)nm, v, s, b, 1U, obs, (double)e, GMT_SVD_EIGEN_NUMBER_CUTOFF);
/* obs (hence alpha) now has the solution for the coefficients based on the first e eigenvalues */
for (k = 0; k < nm; k++)
modnorm += obs[k] * obs[k];
GMT_Report (API, GMT_MSG_INFORMATION, "Used Eigenvalues, Model norm = %d, %g\n", e, modnorm);
Thanks for the nicely maintained and commented source code and easy compilation instructions
Looking at this, could your modnorm value be added as a column to the output file that -E sets? Would be pretty simple for me to add that, possibly via a -E modifier.
The eigenvalues is not what I need, I need the alphas, dependent on n_use.
From C+f I get what I think would be what be the singular values (The diagonal elements of S, from G = USV^T; they are in the 1e18 - 1e0 range. I would expect the alphas to span at most 5 or 6 orders of magnitude). Moreover, I only get one set of e values (e = number of data points = number of model parameters), not e times e.
Yes, that would probably be the easiest one. In fact, two extra columns would be nice: model norm (as modnorm, above) and data norm (with the weights applied just as in the misfit function).
set (CMAKE_INSTALL_PREFIX /home/myusername/.local)
set (GMT_INSTALL_MODULE_LINKS TRUE)
and
$ which greenspline
/home/myusername/.local/bin/greenspline
This previously worked well with my local commits and greenspline was the (almost) only thing recompiling after I pulled in your changes. Am I doing anything incorrectly?
My output is this (note version and synopsis. Not I have nothing to run it on so cannot say if you are wrong but the code looks for +n and strips it off like we do in lots of places in Gmt;
