# Efficiently fill sparse matrix involving BSpline

## Context

In order to implement regularised fitting in 2 or 3 dimensions (as is done say here) using BSplineBasis one needs to evaluate a basis over a set of data.

For the sake of providing a concrete example,

RX = 8; RY = 8;
F[x_, y_] := Sin[0.5 y] Cos[0.9 x];
data = N[Flatten[Table[{x, y, F[x, y]}, {x, -RX, RX, 1}, {y, -RY,RY, 1}], 1]];
data = data + RandomVariate[NormalDistribution[0, 0.1], {Length[data], 3}];
knots = Range[-RX - 2, RX + 2];


A choice of basis using BSplineBasis could be

basis = Flatten@
Table[BSplineBasis[{3, knots}, i, x] *
BSplineBasis[{3, knots}, j, y], {i, 0, 2 RX}, {j, 0, 2 RY}];


Then filling naively the matrix is done as follows:

ff = Function[{x, y}, basis // Evaluate];
a = ff @@ # & /@ (Most /@ data);// AbsoluteTiming


(* 0.245 sec *)

 a // ByteCount


(* 2081056 *)

whereas

  a//SparseArray // ByteCount


(* 68784 *)

## Question

How to make filling the sparse matrix a more efficient in Mathematica so as to be able to deal with large data sets?

I am looking for a solution which minimises both the ByteCount and and AbsoluteTiming.

By construction the matrix a is very sparse.

a // MatrixPlot


so I am guessing it must be possible to fill sparsely the matrix?

PS: The answer would be of interest for order order of BSplines. E.g.

 basis = Flatten@
Table[BSplineBasis[{1, knots}, i, x] *
BSplineBasis[{1, knots}, j, y], {i, 0, 2 RX}, {j, 0, 2 RY}];


PS2: why do a // SparseArray // MatrixPlot and a// MatrixPlot differ visually?

Sparse construction:

With[{xOrder = Ordering[Join[data[[All, 1]], knots]],
yOrder = Ordering[Join[data[[All, 2]], knots]]},
With[{xPar = xOrder[[# + 1 ;; #2 - 1]] & @@@ Partition[Ordering[xOrder, -Length[knots]], 2, 1],
yPar = yOrder[[# + 1 ;; #2 - 1]] & @@@ Partition[Ordering[yOrder, -Length[knots]], 2, 1]},
nonzero = Join @@ Outer[Intersection, Union @@@ Partition[xPar, 4, 1],
Union @@@ Partition[yPar, 4, 1], 1];]]

colIndex = Range[Length[basis]];