Solving the system purely numerically (without any symbolic computations) is a bit faster. This can be done by using the signed incidence matrix of the graph (matrix `B` below).

I have to use the `\[Theta]_?VectorQ` pattern in otder to prevent `NDSolve` from starting any symbolic analysis.

 `ParametricNDSolveValue` provides a neat way to collect everyting into a single function. 

Also, there is no point in solving the ODE until `tmax` if the solution is evaluated only at `evaltime`; so I let just return the state at that time, not the whole `InterpolationFunction`.

    ClearAll[\[Lambda], T];
    B = N[IncidenceMatrix[Graph[UpperTriangularize[A]["NonzeroPositions"],DirectedEdges -> True]]];
    X[\[Theta]_?VectorQ, \[Lambda]_] := N[omegadata] - \[Lambda] B . Sin[\[Theta] . B];
    newa = ParametricNDSolveValue[{
      \[Theta]'[t] == X[\[Theta][t], \[Lambda]], 
      \[Theta][0] == \[Theta]start}, \[Theta][T], 
      {t, 0, T}, {\[Lambda], T}];
    
    aa = a[0.2, evaltime]; // AbsoluteTiming // First
    bb = newa[0.2, evaltime]; // AbsoluteTiming // First
    Max[Abs[aa - bb]/Max[Abs[aa]]]

> 1.65085
> 
> 0.019422
> 
> 4.63889*10^-8

This turns out to be almost 100 times faster. Probably a different ODE integrator is applied, which explains the small difference of the results.