bgrimstad / splinter

Library for multivariate function approximation with splines (B-spline, P-spline, and more) with interfaces to C++, C, Python and MATLAB
Mozilla Public License 2.0
418 stars 115 forks source link

SparseLU solver is inapplicable to non-sqare matrices #106

Closed IuriiSorokin closed 6 years ago

IuriiSorokin commented 6 years ago

When fitting a BSpline to N > 100 data points, the SparseLU solver is used. However, if no regularization is chosen, in compute_control_points(...) the matrices A and B are not square, and SparseLU is inapplicable to non-square matrices. The problem can be reproduced with the code below.

I assume, this is not the expected behavior, otherwise please delete the issue.

#include "SPLINTER/bspline.h"
#include "SPLINTER/data_table.h"
#include <math.h>
#include <random>

using namespace SPLINTER;

int main()
{
    auto random_engine = std::default_random_engine();

    // Function, the data points will follow
    const auto func = []( double x ){ return std::sin(x); };

    // Range
    const double x_min = 0;
    const double x_max = 2*M_PI;

    // Scattering of the points w.r.t. the function
    const double sigma = 1;
    auto scattering_distr = std::normal_distribution<double>(0, sigma);

    // Generate data
    const size_t n_points = 150;
    DataTable data_table;
    auto x_distr = std::uniform_real_distribution<double>(x_min, x_max);
    for( size_t i_point = 0; i_point < n_points; ++i_point )
    {
        const double x = x_distr( random_engine );
        const double y = func(x) + scattering_distr(random_engine);
        data_table.add_sample( x, y );
    }

    // Define and fit spline
    const unsigned int y_dim = 1;
    const std::vector<unsigned int> degrees = {3};
    const std::vector<std::vector<double>> knots = { {x_min, x_min, x_min, x_min, x_max, x_max, x_max } };
    BSpline spline( degrees, knots, y_dim );
    spline.fit( data_table );
}
bgrimstad commented 6 years ago

Hi,

Thank you for reporting and supplying an example.

I failed to catch this requirement of SparseLU last time I updated the _compute_controlpoints function. I have committed a fix in ce953636453941b98e6bd4fdf166dae8068d30f9. We now use the SparseQR solver if the left-hand side matrix is non-square. I have also added your example to test this case (see test/examples/fit_random_1d_spline.cpp). You can run the tests in test/examples by supplying the argument [examples] to splinter-test.

I will proceed with closing this issue.