Efficiently fill sparse matrix involving BSpline

Question

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?

Answer 1

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]];

a2 = SparseArray[Join @@ MapThread[Thread[Thread[{#2, #3}] -> 
          Function[{x, y}, #] @@@ data[[#2, {1, 2}]]] &, {basis, nonzero, colIndex}]];

a == a2

True

The 1st line sorts the x-values along with the knots.
The 3rd line determines where the knots ended up and partitions the data indices according to which 2 knots the x-value is between.
The support of each cubic basis function equals the union of 4 subintervals.