How can I prevent conversion of a SparseArray to a DenseArray?

Question

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.

Answer 1

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