1
$\begingroup$

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

enter image description here

enter image description here

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

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Nov 28, 2023 at 14:09
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Nov 28, 2023 at 14:12
  • $\begingroup$ @MarcoB Thanks, if dq=10, the least squares for sparse matrix “joinjoin‘’ (about 3 seconds)is also slower than normal "joinjoin"(0.283116 seconds) $\endgroup$ Commented Nov 28, 2023 at 14:15

1 Answer 1

5
$\begingroup$

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

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

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.