# Speed up selecting positive eigenvalues repeatedly

I have a smallish (e.g. 2x2 or 4x4, but ideally up to 10x10) non-symmetric square matrix $$\mathbf{A}(x)$$. I need to define a function $$f(x)$$ which is the sum of the eigenvalues with positive real part of $$\mathbf{A}(x)$$, for $$x \in [x_a,x_b]$$.

I then wish to use this function as a normalisation factor within a system of differential equations, $$\mathbf{y}= (\mathbf{Q}(x) - f(x) \cdot \mathbf{I}) \,\mathbf{y}$$.

For a small initial example, for the 2x2 case $$\mathbf{Q}=\mathbf{A}$$.

A = {{0, 1}, {-1 + x^2, 0}};
Q = A;


The eigenvalues here become imaginary for $$|x|<1$$, so I need to be able to detect those transitions cleanly.

Plot[Evaluate@Eigenvalues[A], {x, -4, 4}]


xa = -4; xb = 4;
yVector = {y[1][x], y[2][x]};
ICs = {y[1][xa] == 0, y[2][xa] == 1};

Clear[eee, totalPosEigs];
eee[z_] = Eigenvalues[N[A/.x->z]];
totalPosEigs[z_?NumericQ] :=
totalPosEigs[z] = (Sow[z]; Total@Select[eee[z], Re[#] > 0 &]);
yeqn = Thread[D[yVector, x] == (Q - totalPosEigs[x]*IdentityMatrix[Length[Q]]).yVector];
NDSolve[{yeqn, ICs},
Array[y, {Length[yVector]}], {x, xa, xb}]; // AbsoluteTiming
(* {0.01045, Null} *)


This works, and is fast enough for the 2x2 case. However, as the matrix $$\mathbf{A}$$ gets larger (here 4x4), this approach gets a bit too clunky (and Q is now 6x6):

A = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1},
{-(1/ϵ^4), -((4 Cos[x])/ϵ^2), -((4 Sin[x])/ϵ^2), 0}} /. ϵ -> 1/10;
Q = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0}, {-220 Cos[x], -220 Sin[x], 0, 0, 1, 0},
{0, 0, 0, 0, 1, 0}, {10000, 0, 0, -220 Sin[x], 0, 1}, {0, 10000, 0, 220 Cos[x], 0, 0}};


The eigenvalue structure is also more complicated in this case:

Plot[Evaluate@Re@Eigenvalues[A], {x, xa, xb}]


xa = -4; xb = 4;
yVector = Through[Array[y, {Length[Q]}][x]];
ICs = {y[1][xa] == 0, y[2][xa] == 0, y[3][xa] == 0, y[4][xa] == 0,
y[5][xa] == -1, y[6][xa] == 0};
Clear[eee, totalPosEigs];
eee[z_] = Eigenvalues[N[A /. x -> z]];
totalPosEigs[z_?NumericQ] :=
totalPosEigs[z] = (Sow[z]; Chop@Total@Select[eee[z], Re[#] > 0 &]);
yeqn = Thread[D[yVector, x] == (Q - totalPosEigs[x]*IdentityMatrix[Length[Q]]).yVector];
AbsoluteTiming[{sol, {pts}} =
Reap[NDSolve[{yeqn, ICs}, Array[y, {Length[yVector]}], {x, xa, xb}]];]
{0.229799, Null}


Here totalPosEigs has been evaluated 567 times (even with caching), each evaluation is fast but they are adding up to contribute most of the time of the calculation. If I was to just use the sum of the absolute value of the eigenvalues for instance it takes a quarter of the time, so the Select is a large part of the slowdown.

Its worth noting that the change points may have an imaginary part when they become positive, and the matrix may contain interpolation functions of $$x$$.

• To clarify: do you need the sum of the positive eigenvalues (as your text says: positive real part and zero imaginary part) or the sum of the eigenvalues with positive real part (as your code says: positive real part and arbitrary imaginary part)? – Roman May 15 at 9:09
• @Roman, sorry, positive real part yes. I'll clarify that, thanks. – KraZug May 15 at 10:17
• There is no need to Sow/Reap the points at which totalPosEigs is called, as these are stored as DownValues anyway. You can extract them with pts = Cases[DownValues[totalPosEigs], RuleDelayed[_[totalPosEigs[x_]], y_] /; NumericQ[x] && NumericQ[y] :> {x, y}] after the code finishes. – Roman May 15 at 12:17

Hmmm. I would have set up this system as follows and this way, it is about 5 times as fast:

A = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {-(1/ϵ^4), -((4 Cos[x])/ϵ^2), -((4 Sin[x])/ϵ^2), 0}} /. ϵ -> 1/10;
Q = {{0, 1, 0, 0, 0, 0}, {0, 0, 1, 1, 0, 0}, {-220 Cos[x], -220 Sin[x], 0, 0, 1, 0}, {0, 0, 0, 0, 1, 0}, {10000, 0, 0, -220 Sin[x], 0, 1}, {0, 10000, 0, 220 Cos[x], 0, 0}};

n = Length[Q];
cA = With[{code = N@A}, Compile[{{x, _Real}}, code, CompilationTarget -> "C"]];
cQ = With[{code = N@Q}, Compile[{{x, _Real}}, code, CompilationTarget -> "C"]];

sysmat[x_?NumericQ] := With[{
λ = (1. - UnitStep[-Re[#]]).# &@Eigenvalues[cA[x]]
},
cQ[x] - DiagonalMatrix[ConstantArray[λ, n]]
];

xa = -4; xb = 4;
Y0 = {0, 0, 0, 0, -1, 0};
Ysol = NDSolveValue[{
Y'[x] == sysmat[x].Y[x],
Y[xa] == Y0
},
Y, {x, xa, xb}
]; // AbsoluteTiming // First


0.043133

Most essential steps:

• Compiling definitions of A and Q.
• Teplacing Select by UnitStep and Total by Dot.
• Performing all relevant computation in sysmat with scoped variables so that memoization is actually not required.

OP asked whether this would work also with ParametricNDSolveValue if A and Q depend on a further parameter, say z. My first idea was to add z as a further argument sysmat to sysmat, but upon calling ParametricNDSolveValue, this resulted in an error message

"Dependent variables {Y,sysmat[x,z]} cannot depend on parameters {z}."

I hve no idea for what reason this happens. I consider it a bug.

One can circumvent this by submitting also Y[x] as variable to the function F on the right hand side of the ODE:

F[x_?NumericQ, z_?NumericQ, Y_?VectorQ] :=
With[{λ = (1. - UnitStep[-Re[#]]).# &@Eigenvalues[cA[x]]},
(cQ[x] - DiagonalMatrix[ConstantArray[λ, n]]).Y
];

Y0 = {0, 0, 0, 0, -1, 0};
Ysol = ParametricNDSolveValue[
{
Y'[x] == F[x, z, Y[x]],
Y[xa] == Y0
},
Y,
{x, xa, xb},
z
];


(I am well aware that the parameter z does not effect the result.)

• Thank you Henrik. I don't think the compilation is making any real difference here, but the UnitStep/Dot magic instead of Select makes such a difference. – KraZug May 15 at 12:54
• Actually, compiling the functions cA and cQ here takes ~0.1 seconds each, the compilation is helping the speed (when there are no interpolation functions in the matrix, else compiling needs MainEvaluate anyway) but at the trade-off of the time to compile it. I can probably rewrite my code to only need to compile once for multiple parameters. – KraZug May 15 at 15:05
• Yes, of course: You may add just further arguments to the compiled functions to supply further parameters. – Henrik Schumacher May 15 at 15:16
• Thanks Henrik, i rewrote my code to only compile once for all x and another parameter and to use the vectorized form. – KraZug May 16 at 18:59
• You're welcome! – Henrik Schumacher May 16 at 19:00