6
$\begingroup$

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

Mathematica graphics

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?

$\endgroup$
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ It seems at least 10X faster thanks $\endgroup$ – chris Jul 2 '18 at 18:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.