Also as far as I'm aware dynamic memory allocation in C# is fast (I read it's just a pointer increment, not sure how it works).

Dynamic Memory Allocation is fast in .Net 4+ and in newer versions of Mono. Unfortunately Unity still uses Mono 2.6 and that has a pretty bad garbage collector, so you have to be careful of memory allocation and especially destroying if you need good performance. I recommend you read the following article series: http://www.gamasutra.com/blogs/WendelinReich/20131109/203841/C_Memory_Management_for_Unity_Developers_part_1_of_3.php

But I still think your best solution would be to slice up your algorithm using Coroutines. Heavy calculations are best spread over multiple frames, if you don't need precise real time adjustments. You don't need to calculate the positions of every agent 30+ times per second.

If you don't know about Unitys Coroutines, I gladly provide a few articles about it.

When using the following steering behaviors I had the results below - using 500 boids.
Time to build K-D tree: 2.34ms. So this is not the problem.

Time to update boid steerings (all 500 boids, average time of N frames)
Behaviors used: Cohesion, Separation, VelocityMatch, ObstacleAvoidance => 46.77ms
Behaviors used: Cohesion, Separation, VelocityMatch, ObstacleAvoidance => 39.55ms
Behaviors used: Cohesion, Separation, VelocityMatch, ObstacleAvoidance => 24.44ms

2.36ms is 14% of a 60Hz frame -- quite a large percentage of time! If every task in your game took this long, you could only have 7 different tasks making up your entire game logic (and still run at 60fps).

A good way to get some perspective when profiling these numbers, is to pick a framerate that you'd ideally like to have your game running at (e.g. 60fps or 30fps) and then convert the frame-times (ms) into percentages (e.g. a 60Hz frame gives you a budget of 16.667ms per frame):

Task                              Time      60Hz%  30Hz%
K-D Tree                       => 2.34ms  / 14%  / 7%
VelocityMatch                  => 7.22ms  / 43%  / 22%
Cohesion                       => 22.33ms / 134% / 67%
Separation + ObstacleAvoidance => 17.22ms / 103% / 52%
                                  49.11ms / 294% / 148% < totals

I wouldn't be surprised if it weren't possible to reimplement your entire boids system to process those 500 objects in under 2ms total

Hodge is quite correct, that's a pretty horrible update per frame to be hitting. I also don't see a performance time for:

int neighbors = (maxNeighbors == 0 ? 0 : alglib.kdtreequeryaknn(kdt, point, maxNeighbors, aknnApproxVal));
alglib.kdtreequeryresultsxy(kdt, ref results);

Tree searches are expensive, even via kd-tree's. The first step of that search is going down the tree looking for the closest tree node, then it has to search up and out for the neighbors. If nothing else, this is blowing huge chunks of memory bandwidth and most likely showing up as most of your per boid update. Doing this from within every boid update, almost guarantee'd to be O(log n) or worse in cost.

My suggestion for this would be to move the creation of the list of neighbors out to a system designed to linearize the cost of maintaining such information. You can generally google for: "Spatial Awareness" with things like "c++ library" and get some starting points. I wrote a library which I tested with flocking a while back and even did a vid capture:

Anyway, you are highly unlikely to be able to optimize a per-object solution such as this beyond a certain point. If you follow all the current suggestions, you might drop 5-10ms off the frame times but you'll still have all the same problems, just pushed out a bit further. Try some Googles for "AI Spatial Awareness" including "c++ library" and other trailers and you might find some ideas. You'll be looking for items which supply "aware of" lists. They can be implemented in many ways though I favor sweep and prune due to it having the fewest edge cases and the most consistent per frame processing time no matter what is going on.

It's pretty standard to not run expensive stuff every frame, and AI is often expensive. =)

About co-routines: Seems I do not have access to startCoroutine because the class that handles the updates does not derive from MonoBehavior. However I did it myself: I manually update only n boids every frame (basically what I do is set a flag in each behavior that indicates whether it will use its old value or calculate a new one).

Well, you can start a coroutine from another class (that is a monobehaviour). It just needs to return IEnumerator, that's all. There are some tricks to start coroutines from non-monobehaviour classes as well, but that is pretty deep stuff and I don't think it will help at all in your situation. If you have found a way to spread your calculation over multiple frames - good.

Another thing: If I set more boids (5000 for example ) I still get very long update times even when I still update 100 boids per frame. This is because those that get updaded still do queries on a huge k-d tree... Can I do something about this? What I was thinking was having something that "remembers" where the closest neighbors were so that next time they won't be far away. But how can I know this without re-calculating everything?

And well, for that you can use a heuristic algorithm that uses the last known position, and velocity vector of an agent and tries to "guess" where it will be after a few frames. Of course, the longer you "guess" the more errors in calculation you will have. So you need to balance your update times well, so that you still got good flocking behaviour.

Whilst everyone is correct in terms of addressing memory allocations and re-engineering your overall algorithms being used will give the best improvements, I thought it odd that no-one had pointed out that you are using a standard distance check instead of a squared distance check. I.e. instead of using Vector3.Distance() use Vector3.sqrMagnitude() and provide your MaxDistance as a squared value. Though these days i've no idea of the relative performance benefits of doing this, probably not as much as it used to, but still worth doing IMHO.

I'd also maybe look into where you convert the KD Tree results into Unity types (Vector3). For example you convert vel but its only used if several conditional checks pass, would it not be more prudent to do the conversion only when you needed to add it to the avgVeL?

Indeed I would maybe question why you convert to Unity types at all within the function, leave the vectors as doubles (unless there is a tangible benefit in casting to and working with floats) and only convert to Unity Vector3 the avgCenter and avgVeL at the end of the function or in the calling function. It would mean writing your own sqrDistance and dotproduct methods, but they are easy enough, just don't forget to convert your character.pos from Vector3 to a double array in the calling function. In other words its probably far better to keep all your data in one format/type and convert at the end, instead of converting on the fly, especially when some of the data may not even be used.

However as I said at the beginning you probably find better gains through changing your approach or algorithms, but the few optimisations above are pretty straightforward and simple to implement.

For repeating an action infrequently like updating your KDtree you could use InvokeRepeating() in Unity.

Its not clear to me but are you calling GetNeighborhoodInfo() for each boid for each Cohesion, Separation, VelocityMatch? If so shouldn't you simply be sharing the results between all three behaviours and just calculate it once per AI Tick using the highest K value you need?

Finally a few points about Unity. If you have a spare email address you should be able to sign up for a 30 day demo of Pro, getting you access to the profiler. Alternatively you could try emailing Unity Technologies and asking for a Pro demo key, explaining that you are using it for your thesis and supply evidence, as they are quite a nice company and might be willing to help you out - at least it doesn't hurt to try.