# Dot Product of Matrices in NDSolve

How to efficiently solve a system of differential equations involving matrix dot product?

The example below is solved extremely fast when the array structure is flattened out. Is there a way to reach similar speed preserving the array structure?

x = 4; y = 2;
(ini = SparseArray[{1, 1, 1} -> 1., {x, 1, y}]) // MatrixForm
(a = SparseArray[
Flatten[Table[{x1, x2, y1, y2} -> N[ x1/(x2 y1), 2], {y2, y}, {y1,
y}, {x2, x}, {x1, x}], 4], {x, x, y, y}]) // MatrixForm  Flatten

AbsoluteTiming[sol1 = NDSolve[
{
f'[t] == Normal[ArrayFlatten[a].Flatten[f[t]]],
f == Normal[Flatten[ini]]
}
,f,
{t, 0, 2}
];]
(*{0.007967, Null}*)


No Flatten

dodo[m1_, m2_, pos_] := TensorTranspose[Activate[TensorContract[Inactive[TensorProduct][m1,m2], pos]], {1,3, 2}];
AbsoluteTiming[sol2 = NDSolve[
{
f'[t] == Normal[dodo[a, f[t], {{2, 5}, {4, 7}}]],
f == Normal[ini]
}
,f,
{t, 0, 2}
];]
(*{0.195254, Null}*)


Both solutions are equal

Round[Flatten[f /. sol1], 10^-5] == Round[Flatten[f /. sol2], 10^-5]
(*True*)

• Why couldn't you just re-partition the flattened solution afterwards? Nov 2 '17 at 15:56
• The original problem is bigger and easier to read in the array form. Of course, re-partitioning would be plan B. Nov 2 '17 at 16:12
• "The original problem is bigger and easier to read in the array form." - I don't doubt that. I'm just saying that this seems to be a situation where you have to choose between readability and efficiency. Nov 2 '17 at 16:16

You can improve timings by giving dodo an argument restriction, so that it doesn't Activate too early:

Clear[dodo];
dodo[m1_, m2_List, pos_] := TensorTranspose[
Activate[TensorContract[Inactive[TensorProduct][m1,m2],pos]],
{1,3,2}
];
AbsoluteTiming[
sol2 = NDSolve[
{f'[t]==Normal[dodo[a,f[t],{{2,5},{4,7}}]],f==Normal[ini]},
f,
{t,0,2}
];
]


{0.029865, Null}

Clear[dodo];
dodo[m1_, m2_, pos_] := TensorTranspose[
Activate[TensorContract[Inactive[TensorProduct][m1,m2],pos]],
{1,3,2}
];
AbsoluteTiming[
sol3 = NDSolve[
{f'[t]==Normal[dodo1[a,f[t],{{2,5},{4,7}}]],f==Normal[ini]},
f,
{t,0,2}
];
]

Round[Flatten[f/.sol3],10^-5]==Round[Flatten[f/.sol2],10^-5]


{0.223253, Null}

True

This is still quite a bit slower than working with flattened objects, though. It is possible to improve things further, but not enough to be competitive with flattened objects.

• Thanks! That helps but with higher dimensions things get slow again.. I guess there is no way around Flatten. Just out of curiosity, how would you improve the array version? Nov 2 '17 at 17:56