3
$\begingroup$

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

enter image description here

enter image description here

Flatten

AbsoluteTiming[sol1 = NDSolve[
{
 f'[t] == Normal[ArrayFlatten[a].Flatten[f[t]]],
 f[0] == 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[0] == Normal[ini]
}
,f,
{t, 0, 2}
];]
(*{0.195254, Null}*)

Both solutions are equal

Round[Flatten[f[1] /. sol1], 10^-5] == Round[Flatten[f[1] /. sol2], 10^-5]
(*True*)
$\endgroup$
  • 1
    $\begingroup$ Why couldn't you just re-partition the flattened solution afterwards? $\endgroup$ – J. M. is away Nov 2 '17 at 15:56
  • $\begingroup$ The original problem is bigger and easier to read in the array form. Of course, re-partitioning would be plan B. $\endgroup$ – tukan Nov 2 '17 at 16:12
  • 1
    $\begingroup$ "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. $\endgroup$ – J. M. is away Nov 2 '17 at 16:16
3
$\begingroup$

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[0]==Normal[ini]},
        f,
        {t,0,2}
    ];
]

{0.029865, Null}

Compare to your version:

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[0]==Normal[ini]},
        f,
        {t,0,2}
    ];
]

Round[Flatten[f[1]/.sol3],10^-5]==Round[Flatten[f[1]/.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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – tukan Nov 2 '17 at 17:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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