Speed up selecting positive eigenvalues repeatedly

Question

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$.

Answer 1

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.

Addendum

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.)