As of the time of writing, the current version essentially:
-
Applies the 1-D knot construction algorithm to each dimension
The input data for the independent variable is rectangular, so this only needs to be done once per dimension.
-
Applies the 1-D coefficient solution algorithm to each line in each dimension
The algorithm is based upon solving a linear system of equations, Ax=b, where A depends only on the independent variable sample points (and the constructed knots) and b depends on the value at the point. For the ith dimension of a general N-dimensional spline, there are different observations for each independent variables. For example, in the two dimensional case, using x and y for the independent variables, for same set of samples at x_1, x_2, ..., x_n, there are is an observation for each y_1, y_2, ..., y_m. It turns out we can apply the same algorithm to every line of the tensors of sampled data to compute the coefficients for the N-dimensional spline.
As mentioned in #27, this approach may not be the best. At some point it would nice to do some numerical analysis and/or experimentation to compare the current algorithm to a new approach(es). Two possibilities are:
-
Use the same general approach as above, but solve for all the (parallel) lines of each dimension simultaneously. This should reduce computation time by factoring the constant A only once. The SciPy impelemtnation picks out a particular LAPACK solver (gbsv) to take advantage of the banded structure. It looks like it can handle multiple right-hand side b vectors.
-
Solve for all coefficients simultaneously, similar to the ndspline make_lsq_spline approach. This has the possibility of changing the numerical complexity, although I would need to look into how to analytically compare multiple sequential solutions compared to the single solution. Experimentally, we can use timing and we can compare to SciPy's implementations for 1-D and 2-D cases for coefficient accuracy (I think this is equivalent to numerical stability).
As of the time of writing, the current version essentially:
Applies the 1-D knot construction algorithm to each dimension
The input data for the independent variable is rectangular, so this only needs to be done once per dimension.
Applies the 1-D coefficient solution algorithm to each line in each dimension
The algorithm is based upon solving a linear system of equations,
Ax=b, whereAdepends only on the independent variable sample points (and the constructed knots) andbdepends on the value at the point. For theith dimension of a general N-dimensional spline, there are different observations for each independent variables. For example, in the two dimensional case, usingxandyfor the independent variables, for same set of samples atx_1, x_2, ..., x_n, there are is an observation for eachy_1, y_2, ..., y_m. It turns out we can apply the same algorithm to every line of the tensors of sampled data to compute the coefficients for the N-dimensional spline.As mentioned in #27, this approach may not be the best. At some point it would nice to do some numerical analysis and/or experimentation to compare the current algorithm to a new approach(es). Two possibilities are:
Use the same general approach as above, but solve for all the (parallel) lines of each dimension simultaneously. This should reduce computation time by factoring the constant
Aonly once. The SciPy impelemtnation picks out a particular LAPACK solver (gbsv) to take advantage of the banded structure. It looks like it can handle multiple right-hand sidebvectors.Solve for all coefficients simultaneously, similar to the ndspline
make_lsq_splineapproach. This has the possibility of changing the numerical complexity, although I would need to look into how to analytically compare multiple sequential solutions compared to the single solution. Experimentally, we can use timing and we can compare to SciPy's implementations for 1-D and 2-D cases for coefficient accuracy (I think this is equivalent to numerical stability).