# Why does using a sparse matrix for least squares instead slow down the speed?

Here is the code. The matrix appears to be relatively sparse in terms of its visual representation, but why does using least squares result in significantly slower speed compared to the previously dense matrix? (The code appears to be quite lengthy, but it only takes about two minutes to reproduce. I currently do not have a way to reduce it to a minimal example.) May I ask why it causes the slowdown, and if there is any solution to it? Thanks, Thanks.

test[t_, x_, y_] :=
Sin[2 \[Pi] x] Sin[
2 \[Pi] y] (2 Cos[Sqrt[(2 \[Pi]/5)^2 + (2 \[Pi]/4)^2] t 10] +
Sin[Sqrt[(2 \[Pi]/5)^2 + (2 \[Pi]/4)^2] t 10]);
deltau[t_, x_, y_] =
0.01 D[D[test[t, x, y], t], t] - 1/25 D[D[test[t, x, y], x], x] -
1/16 D[D[test[t, x, y], y], y];
deltaT[t_, x_, y_] = D[test[t, x, y], t];
deltau =
With[{code = deltau[t, x, y]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]];
deltaT =
With[{code = deltaT[t, x, y]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]];
test = With[{code = test[t, x, y]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]];

(*CloseKernels[];
LaunchKernels[6];*)
(*\$Pre=Function[Null,MemoryConstrained[#,1000000000-MemoryInUse[]],\
HoldAll];*)

dimension = 3;
n = 2.;
dn = (n)^dimension;(*和维数有关*)
dj = 50.;
(*dq=5.;*)
collocation = Tuples[Table[(i)/(n), {i, 0, n - 1}], dimension];
random1 = RandomReal[{-1, 1}, {dn, dj}];
random2 = RandomReal[{-1, 1}, {dn, dj}];
random3 = RandomReal[{-1, 1}, {dn, dj}];
random4 = RandomReal[{-1, 1}, {dn, dj}];
random5 = RandomReal[{-1, 1}, {dn, dj}];

psi[i_, t_, x_, y_] := 1.(*f[i,1,x]*f[i,2,y]*);
li[i_, t_, x_,
y_] := {(t - (collocation[[i]][[1]] + 1/(2 n)))/
Abs[1/(2 n)], (x - (collocation[[i]][[2]] + 1/(2 n)))/
Abs[1/(2 n)], (y - (collocation[[i]][[3]] + 1/(2 n)))/
Abs[1/(2 n)]};
phi[i_, j_, t_, x_, y_] :=
Tanh[{random1[[i, j]], random2[[i, j]], random3[[i, j]]} .
li[i, t, x, y] + random5[[i, j]]];
df[t_, x_, y_] :=
Table[phi[i, j, t, x, y]*psi[i, t, x, y], {i, 1, dn}, {j, 1, dj}] //
Flatten;
ddt[t_, x_, y_] = D[D[df[t, x, y], t], t]; // Timing
ddx[t_, x_, y_] = D[D[df[t, x, y], x], x];
ddy[t_, x_, y_] = D[D[df[t, x, y], y], y];

dd3[t_, x_, y_] =
1/100 ddt[t, x, y] - 1/25 ddx[t, x, y] - 1/16 ddy[t, x, y];

xdim = IntegerPart[dn dj];

funComptotal =
With[{code = N[dd3[t, x, y]]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]]; // AbsoluteTiming
funCompdf =
With[{code = N[df[t, x, y]]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]]; // AbsoluteTiming
funCompDt =
With[{code = N[D[df[t, x, y], t]]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]];
funCompDx =
With[{code = N[D[df[t, x, y], x]]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]];
funCompDy =
With[{code = N[D[df[t, x, y], y]]},
Compile[{{t, _Real}, {x, _Real}, {y, _Real}}, code,
CompilationTarget -> "C"]];

dq = 20;
dQ = dq + 1;
omegapoint =
Table[Tuples[{Table[i/(n dq), {i, (j - 1) dq, j dq}],
Table[i/(n dq), {i, (k - 1) dq, k dq}],
Table[i/(n dq), {i, (p - 1) dq, p dq}]}], {j, 1, n}, {k, 1,
n}, {p, 1, n}];
omegapointbd5 =
Table[Tuples[{Table[(i - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(j - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(k - 1)/n, {z, 1}]}], {i, n}, {j, n}, {k, n}];(*向t-x投影0*)
(*omegapointbd1=Table[{collocation[[i]][[1]]+(q-1)/(n \
dq),0.},{i,1,dn},{q,1,dQ}];*)
omegapointbd6 =
Table[Tuples[{Table[(i - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(j - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(k)/n, {z, 1}]}], {i, n}, {j, n}, {k, n}];(*向t-x投影1*)
(*omegapointbd2=Table[{collocation[[i]][[1]]+(q-1)/(n \
dq),1.},{i,1,dn},{q,1,dQ}];*)
omegapointbd3 =
Table[Tuples[{Table[(i - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(j - 1)/n, {z, 1}],
Table[(k - 1)/n + (z - 1)/(n dq), {z, dQ}]}], {i, n}, {j, n}, {k,
n}];(*向t-y投影0*)

omegapointbd4 =
Table[Tuples[{Table[(i - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(j)/n, {z, 1}],
Table[(k - 1)/n + (z - 1)/(n dq), {z, dQ}]}], {i, n}, {j, n}, {k,
n}];(*向t-y投影1*)

omegapointbd1 =
Table[Tuples[{Table[(i - 1)/n, {z, 1}],
Table[(j - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(k - 1)/n + (z - 1)/(n dq), {z, dQ}]}], {i, n}, {j, n}, {k,
n}];(*向x-y投影0*)
omegapointbd2 =
Table[Tuples[{Table[(i)/n, {z, 1}],
Table[(j - 1)/n + (z - 1)/(n dq), {z, dQ}],
Table[(k - 1)/n + (z - 1)/(n dq), {z, dQ}]}], {i, n}, {j, n}, {k,
n}];(*向x-y投影1*)

tabledfpoint = SparseArray[Table[Table[0., {xdim}], {(n dQ)^3}]];
tabledfpointbd1 = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
tabledfpointbd2 = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
tabledfpointbd3 = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
tabledfpointbd4 = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
tabledfpointbd5 = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
tabledfpointbd6 = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];

Do[tabledfpoint[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^3 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^3, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funComptotal @@@ omegapoint[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^3 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^3, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]], {i, 1, n}, {j, 1,
n}, {k, 1, n}]; // AbsoluteTiming
Do[If[ k == 1,
tabledfpointbd1[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funCompdf @@@ omegapointbd5[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]]], {i, 1, n}, {j, 1,
n}, {k, 1, n}];
Do[If[ k == n,
tabledfpointbd2[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funCompdf @@@ omegapointbd6[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]]], {i, 1, n}, {j, 1,
n}, {k, 1, n}];

Do[If[ i == 1,
tabledfpointbd5[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funCompdf @@@ omegapointbd1[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]]], {i, 1, n}, {j, 1,
n}, {k, 1, n}];
(*Do[If[ i==n,tabledfpointbd4[[((i-1)n+(j-1))dQ+1;;((i-1)n+j) \
dQ,((i-1)n+(j-1)) dj+1;;((i-1)n+j) \
dj]]=(Join@@Join@@Table[funCompDt@@@omegapointbd3[[i,j]],{i,1,n},{j,1,\
n}])[[((i-1)n+(j-1))dQ+1;;((i-1)n+j) dQ,((i-1)n+(j-1)) \
dj+1;;((i-1)n+j) dj]]],{i,1,n},{j,1,n}];*)
Do[If[ i == 1,
tabledfpointbd6[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funCompDt @@@ omegapointbd1[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]]], {i, 1, n}, {j, 1,
n}, {k, 1, n}];

Do[If[ j == 1,
tabledfpointbd3[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funCompdf @@@ omegapointbd3[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]]], {i, 1, n}, {j, 1,
n}, {k, 1, n}];
Do[If[ j == n,
tabledfpointbd4[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]] = (Join @@
Join @@ Join @@
Table[funCompdf @@@ omegapointbd4[[i, j, k]], {i, 1, n}, {j,
1, n}, {k, 1,
n}])[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dj]]], {i, 1, n}, {j, 1,
n}, {k, 1, n}];

table0t = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
table1t = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
table0x = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
table1x = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
table0y = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];
table1y = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]];

Do[If[i !=
1, {table0t[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 2) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 2) n^2 + (j - 1) n + (k)) dj]] = -(funCompdf @@@
omegapointbd1[[i, j,
k]])[[;; , ((i - 2) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 2) n^2 + (j - 1) n + (k)) dj]],
table0t[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]] = (funCompdf @@@
omegapointbd1[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]],
table1t[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 2) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 2) n^2 + (j - 1) n + (k)) dj]] = -(funCompDt @@@
omegapointbd1[[i, j,
k]])[[;; , ((i - 2) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 2) n^2 + (j - 1) n + (k)) dj]],
table1t[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]] = (funCompDt @@@
omegapointbd1[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]]}], {i, n}, {j,
n}, {k, n}];
Do[If[j !=
1, {table0x[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 2) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 2) n + (k)) dj]] = -(funCompdf @@@
omegapointbd3[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 2) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 2) n + (k)) dj]],
table0x[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]] = (funCompdf @@@
omegapointbd3[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]],
table1x[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 2) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 2) n + (k)) dj]] = -(funCompDx @@@
omegapointbd3[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 2) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 2) n + (k)) dj]],
table1x[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]] = (funCompDx @@@
omegapointbd3[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]]}], {i, n}, {j,
n}, {k, n}];
Do[If[k !=
1, {table0y[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 2)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k - 1)) dj]] = -(funCompdf @@@
omegapointbd5[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 2)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k - 1)) dj]],
table0y[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]] = (funCompdf @@@
omegapointbd5[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]],
table1y[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 2)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k - 1)) dj]] = -(funCompDy @@@
omegapointbd5[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 2)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k - 1)) dj]],
table1y[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n +
k) dQ^2, ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]] = (funCompDy @@@
omegapointbd5[[i, j,
k]])[[;; , ((i - 1) n^2 + (j - 1) n + (k - 1)) dj +
1 ;; ((i - 1) n^2 + (j - 1) n + (k)) dj]]}], {i, n}, {j,
n}, {k, n}];

joinjoin =
Join[tabledfpoint, tabledfpointbd1, tabledfpointbd2,
tabledfpointbd5, tabledfpointbd6, tabledfpointbd3, tabledfpointbd4,
table0t, table0x, table0y, table1t, table1x, table1y];
lambda = N[100/Replace[Max /@ Abs[joinjoin], 0. -> Infinity, {1}]];
joinjoin = lambda*joinjoin;
point = SparseArray[Table[0., First[Dimensions[tabledfpoint]]]];
pointbd1 =
SparseArray[Table[0., First[Dimensions[tabledfpointbd1]]]];
pointbd2 =
SparseArray[Table[0., First[Dimensions[tabledfpointbd2]]]];
pointbd3 =
SparseArray[Table[0., First[Dimensions[tabledfpointbd3]]]];
pointbd4 =
SparseArray[Table[0., First[Dimensions[tabledfpointbd4]]]];
pointbd5 =
SparseArray[Table[0., First[Dimensions[tabledfpointbd5]]]];
pointbd6 = SparseArray[Table[0., First[Dimensions[tabledfpointbd6]]]];

t0t = SparseArray[Table[0., First[Dimensions[table0t]]]];
t1t = SparseArray[Table[0., First[Dimensions[table1t]]]];
t0x = SparseArray[Table[0., First[Dimensions[table0x]]]];
t1x = SparseArray[Table[0., First[Dimensions[table1x]]]];
t0y = SparseArray[Table[0., First[Dimensions[table0y]]]];
t1y = SparseArray[Table[0., First[Dimensions[table1y]]]];

Do[point[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^3 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^3]] =
deltau @@@ omegapoint[[i, j, k]], {i, n}, {j, n}, {k, n}];
Do[If[ k == 1,
pointbd1[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^2]] =
test @@@ omegapointbd5[[i, j, k]]], {i, 1, n}, {j, 1, n}, {k, 1,
n}];
Do[If[ k == n,
pointbd2[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^2]] =
test @@@ omegapointbd6[[i, j, k]]], {i, 1, n}, {j, 1, n}, {k, 1,
n}];
Do[If[ j == 1,
pointbd3[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^2]] =
test @@@ omegapointbd3[[i, j, k]]], {i, 1, n}, {j, 1, n}, {k, 1,
n}];
Do[If[ j == n,
pointbd4[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^2]] =
test @@@ omegapointbd4[[i, j, k]]], {i, 1, n}, {j, 1, n}, {k, 1,
n}];
Do[If[ i == 1,
pointbd5[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^2]] =
test @@@ omegapointbd1[[i, j, k]]], {i, 1, n}, {j, 1, n}, {k, 1,
n}];
Do[If[ i == 1,
pointbd6[[((i - 1) n^2 + (j - 1) n + (k - 1)) dQ^2 +
1 ;; ((i - 1) n^2 + (j - 1) n + k) dQ^2]] =
deltaT @@@ omegapointbd1[[i, j, k]]], {i, 1, n}, {j, 1, n}, {k, 1,
n}];
f = Join[point, pointbd1, pointbd2, pointbd5, pointbd6, pointbd3,
pointbd4, t0t, t0x, t0y, t1t, t1x, t1y]*lambda;

LeastSquares[joinjoin, f]; // AbsoluteTiming
LeastSquares[Normal[joinjoin], Normal[f]]; // AbsoluteTiming


May I ask why it causes the slowdown, and if there is any solution to it? Thanks, Thanks.

• Try generating a smaller set of matrices with similar structure annd random content, and see if you can reproduce the problem, then post that. Please do some troubleshooting on your own before dumping pages of code here. Nov 28, 2023 at 14:09
• You do table0t = SparseArray[Table[Table[0., {xdim}], {n^3 dQ^2}]]; and then fill it with entries in a loop with table0t [[i,j]] = [...]This is not the way to create large sparse matrices. Please read the documentation in order to learn better ways of creating SparseArrays. Nov 28, 2023 at 14:12
• @MarcoB Thanks, if dq=10, the least squares for sparse matrix “joinjoin‘’ (about 3 seconds)is also slower than normal "joinjoin"(0.283116 seconds) Nov 28, 2023 at 14:15

Have a look at the matrix joinjoin: If you look at its the SparseArray box and then on the "+" sign, it tells you that its density is 0.114. That's a very, very high density for many sparse matrix method. I am not sure, what exactly LeastSquares does (it comes with several method options), but I expect that direct methods, e.g., based on sparse QR decomposition (I think that's what is behind Method-> "LSQR" are used), and I am not surprised that they are not efficient here. Surprising to me is, that Method -> "Direct" actually impoves over the default method; but it is still slower than the dense matrix method.

Since you have many fewer variables than equations, I expect that your matrix A = joinjoin is injective. In that case you can solve the normal equation $$A^T A x = A^T f$$ with a dense matrix method instead:

A = joinjoin;
AT = Transpose[A];
x = LinearSolve[Normal[AT.A], AT.f, Method -> "Cholesky"]; // AbsoluteTiming //First


We use a dense matrix method here because Transpose[A].A has dimensions {400,400}.

0.006359

• Thank you very much, using your method has indeed made the speed much faster. Nov 29, 2023 at 12:07
• You're welcome! Nov 29, 2023 at 12:34