8
$\begingroup$

The Problem

I'm trying to implement an efficient Whitker-Eliers smoother in Mathematica. In Matlab, this is a few lines of code (taken from the SI of the above paper):

m = length(data);
E = speye(m); % Sparse identity matrix
D = diff(E,d); % Equivalent to Differences[E,d]
C = chol(E + lambda * D' * D); % lambda is a constant
smoothed_data = C \ (C' \ y); % \ is equivalent to LinearSolve

When I try a similar setup in Mathematica, the sparse identity matrix almost immediately gets converted to a dense one. Using small matrices for demonstration:

identity = IdentityMatrix[5,SparseArray];
diff = Differences[identity,1]; (* This becomes a dense array *)

I've verified that this is the case at every step of the function chain. I've also verified that Matlab keeps everything as sparse arrays throughout the computation, which may be the reason for the significant speed difference I observe, especially as the size of the matrices grow.

Is there a way to tell Mathematica that I want to keep the SparseArray format when I perform a computation? Moreover, is there a performance benefit to doing so?


Edit:

I realized that I failed to include my own code, and performance benchmarks. In Mathematica:

whitakerSmooth[data_, lambda_, order_]:=
  Module[
    {rawData = data, lambda0 = lambda, d = order},
    identity = IdentityMatrix[Dimensions[rawData][[1]],SparseArray];
    diffMatrix = Differences[identity,d];
    product = identity + lambda0 * Dot[diffMatrix//Transpose,diffMatrix];
    chole = CholeskyDecomposition[product];
    smoothedData = LinearSolve[chole,LinearSolve[chole//Transpose,rawData]];
    Return[smoothedData];
  ]

Benchmarking using Timing:

data = RandomReal[{0,1},1000];
Timing[whitakerSmooth[data,10,1]][[1]]
(* 10.9844 seconds *)

Benchmarking in Matlab:

tic; whitsmtest(data,10,1); toc
% Elapsed time is 0.000863 seconds

Edit 2

Digging a little deeper, it appears that the primary issue is the CholeskyDecomposition step.

AbsoluteTiming[identity = IdentityMatrix[1000,SparseArray]][[1]]
(* 0.000017 *)
AbsoluteTiming[diff = Differences[identity]][[1]]
(* 0.0023012 *)
AbsoluteTiming[Differences[identity,2]][[1]]
(* 0.0035748 *)
AbsoluteTiming[product = identity + 10 * Dot[diff//Transpose, diff]][[1]]
(* 0.0173485 *)
AbsoluteTiming[chole = CholeskyDecomposition[product]][[1]]
(* 11.3763 *)
sparseProduct = SparseArray[product];
AbsoluteTiming[CholeskyDecomposition[sparseProduct]][[1]]
(* 11.2335 *)

@Carl Woll pointed out that the CholeskyDecomposition of an exact matrix is much slower. Numericizing lambda gives the following result:

AbsoluteTiming[product = identity + (N@lambda) * Dot[diff//Transpose,diff]][[1]]
(* 0.0177376 *)
AbsoluteTiming[chole = CholeskyDecomposition[product]][[1]]
(* 0.0194817 *)

Adding this step to the function takes the total timing to the following (including a hacky Rust implementation I wrote):

Mathematica (henrik-schumacher's solution) = 0.76 ms
matlab = 0.863 ms
Mathematica (numericize lambda) = 13.8 ms
python = 96.1 ms
Rust (my poor solution) = 103.98 ms
Mathematica (original solution) = 11 s

Edit 3 Updated above table with more timings

Edit 4 In order to make sure I understood @Henrik Schumacher's solution, I implemented his sparse matrix solution for 2nd order differences:

n = 10;
lambda = 10;
(* Generate a goal matrix *)
identity = IdentityMatrix[n];
diff = Differences[identity,2];
goal = identity + lambda * Dot[diff//Transpose,diff];
goal//MatrixForm
(* Generate our value array *)
vals = ConstantArray[lambda,{5 n - 6}];
(* Manually set a bunch of values *)
vals[[1]] = vals[[-1]] = 1. + lambda;
vals[[{2,4,-2,-4}]] = -2. * lambda;
vals[[10;;-10;;5]] = 1. + 6 * lambda;
vals[[{5,-5}]] = 1. + 5 * lambda;
vals[[6;;-9;;5]] = -4. * lambda;
vals[[9;;-6;;5]] = -4. * lambda;
(* This should return True *)
vals == N@DeleteCases[Flatten[goal],0]
(* Generate the column indices *)
ci = Drop[Drop[(Join@@Partition[Range[-1,n+2],5,1])[[3;;-3]],{4}],{-4}];
(* Generate the row pointers *)
rp = Accumulate[Join[{0},{3},{4},ConstantArray[5,n-4],{4},{3}]];
(* Make some toy data *)
rawData = RandomReal[{0,1},n];
(* Make the actual sparse array using the undocumented method *)
array = With[{rawData = {Automatic, {n,n},0.,{1,{rp,Partition[ci,1]},vals}}},SparseArray@@rawData];
(* This should also return true *)
array == N@goal
(* Solve the system *)
S = LinearSolve[array,Method->"Banded"];
smoothed = S[rawData];
(* Plot the result *)
ListPlot[{rawData,smoothed},Joined->True]

This appears to work and is equivalently fast.

$\endgroup$
5
  • 5
    $\begingroup$ You could do: SparseArray[{Band[{1,1}]->-1, Band[{1,2}]->1}, Dimensions[identity]-{1,0}] . identity instead of using Differences. $\endgroup$
    – Carl Woll
    Commented Jul 26, 2023 at 13:17
  • $\begingroup$ This does meaningfully speed up the computation. For an identity matrix with 1000 entries, AbsoluteTiming gives 0.0015234 for your solution vs 0.0026268 for using Differences. This does get much more complicated if I choose to use Differences order 2 though. I believe this is the correct modification of the code: SparseArray[{Band[{1,2}]->-2,Band[{1,1}]->1,Band[{1,3}]->1},Dimensions[identity]-{2,0}].identity $\endgroup$ Commented Jul 26, 2023 at 14:03
  • 3
    $\begingroup$ The slowest part of the Mathematica code is CholeskyDecomposition of an exact matrix. If you use whitakerSmooth[data, 10., 1] or numericize using lambda0 = N @ lambda, then the Mathematica code will be orders of magnitude faster. $\endgroup$
    – Carl Woll
    Commented Jul 26, 2023 at 14:49
  • $\begingroup$ I was just editing my reply to include that the CholeskyDecomposition was the main issue, but I had no idea how to fix it. This takes the execution time from ~11s to 0.01948s! $\endgroup$ Commented Jul 26, 2023 at 14:54
  • 1
    $\begingroup$ FWIW, Rest@identity - Most@identity is a fast way to compute the (first) differences. Nest[Rest[#] - Most[#] &, identity, d] does $d$th differences, though the speed depends on d. $\endgroup$
    – Michael E2
    Commented Jul 27, 2023 at 0:11

1 Answer 1

9
$\begingroup$

Just for the records: For order = 1, the matrix product is just a sum of the identity matrix plus lambda times the graph Laplacian of a linear graph. The graph Laplacian is also a tiny multiple of the 1D finite difference Laplacian on an equilateral grid of Length[data] points. So what this does is to perform a single backward Euler step of heat flow. No wonder that his works well as a smoother!

I am somewhat shocked to see this method being reinvented in the year 2003. This is pretty close to Tai's method. It's insane how many citations this paper has!

Back to the numerics: As others have already pointed out, it is crucial for performance to use machine precision numbers. The matrix is a sparse, positive-definite, and tridiagonal matrix. So the best method of factorization should be the Thomas algorithm available in LinearSolve with the option Method->"Banded". This is just a Cholesky decomposition; but because the matrix is banded, one can get rid of quite some complicated logic (like determining a fill-inn reducing reordering).

Once you have factorized the matrix, don't throw it away. You can reuse the factorization, for example, to smooth the data multiple times. (It should gain "two derivatives" by pass.)

This is btw. the fasted method I found to build the matrix without resorting to Compile or C++ code:

n = 1000000;

data = RandomReal[{0, 1}, n];
\[Lambda] = 10.;

A = Module[{rp, ci, vals},
     rp = Accumulate[Join[{0}, {2}, ConstantArray[3, n - 2], {2}]];
     ci = (Join @@ Partition[Range[0, n + 1], 3, 1])[[2 ;; -2]];
     vals = ConstantArray[-\[Lambda], {3 n - 2}];
     vals[[1]] = vals[[-1]] = 1. + \[Lambda];
     vals[[4 ;; -4 ;; 3]] = 1. + 2. \[Lambda];
     With[{data = {Automatic, {n, n}, 0., {1, {rp, Partition[ci, 1]}, vals}}}, SparseArray @@ data]]; // AbsoluteTiming // First
S = LinearSolve[A, Method -> "Banded"]; // AbsoluteTiming // First
S[data]; // RepeatedTiming // First

0.185926

0.289245

0.00639384

$\endgroup$
6
  • 2
    $\begingroup$ This reinvention reminds me the story of a recent rediscovery of calculus. You probably refer to the same fact by mentioning the Tai method, or is it something different? $\endgroup$
    – yarchik
    Commented Jul 26, 2023 at 16:40
  • $\begingroup$ Typically, the Whitaker-Eliers smoother is run with differentiation order = 2, so it's likely necessary to do 2 passes, though I'm not sure if that's faster than simply initializing the diff matrix with the second order derivative. I must confess that I'm completely baffled by your code. I'm going to have to break it down line by line. One of the primary advantages of the Whitaker-Eliers smoother (and the reason I use it) is that it's trivial to compute the cross-validation error, which allows for "unattended" smoothing, which is quite useful in large datasets $\endgroup$ Commented Jul 26, 2023 at 17:35
  • $\begingroup$ Your solution came in even faster than the matlab implementation! I've edited the original post to include all of the timings. $\endgroup$ Commented Jul 26, 2023 at 17:59
  • 1
    $\begingroup$ @BesselFunct Well, the matrix A has a very simple sparsity pattern which allows me to write down the data directly into CSR format (rp = row pointers, ci = column indices, val = nonzero values). The remainder is to use an undocumented way to initialize a SparseArray from this data. This skips the costs of assembling the matrix in the first place. A baffling trick, but a cheap one. Nothing too fancy. $\endgroup$ Commented Jul 26, 2023 at 18:40
  • $\begingroup$ @yarchik Exactly that's what I was referring to... ^^ $\endgroup$ Commented Jul 26, 2023 at 18:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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