Extrema of a piecewise sinusoidal function

I consider the function

$$f(\theta)=\sum_{k=1}^n|a_k\cos\theta+b_k\sin\theta|$$ of which I need to find the global maximum (and not necessarily the value of $\theta$ that achieves it).

This function is piecewise sinusoidal, with every piece of the form

$$\cos\theta\sum_{k=1}^n(\pm a_k)+\sin\theta\sum_{k=1}^n(\pm b_k),$$ potentially giving the local maximum

$$\sqrt{\left(\sum_{k=1}^n(\pm a_k)\right)^2+\left(\sum_{k=1}^n(\pm b_k)\right)^2}.$$

The combinations of signs of $a_k\cos\theta+b_k\sin\theta$ delimit intervals of a full period (the period is $\pi$), and a local maximum for a given combination is only relevant if the corresponding $\theta$ value belongs to this interval.

So a possible solution is to build all intervals, solve for the $\theta$'s that achieve a maximum, check if they belong to their interval and compute the values of $f$. As I need to repeat this process million times, what I am after is a shortcut to avoid an $O(n^2)$ computational load ($n$ intervals and cost $O(n)$ of the function evaluation).

Can we find the global maximum more efficiently ?

Sample case:

Update:

If we order the terms by increasing $\theta_k$, where $\theta_k$ is the first root of $a_k\cos\theta+b_k\sin\theta=0$, the update of the sums $\Sigma(\pm a_k)$, $\Sigma(\pm b_k)$ when we cross $\theta_k$ only involves two terms, hence is done in constant time. And the computation of the maximum in an interval is also done in constant time. So we can achieve $O(n)$ after sorting on the angles.

I see a $O(n)$ algorithm defined this way: first write \begin{equation} a_k \cos\theta + b_k\sin\theta = c_k \sin(\theta + \varphi_k) \end{equation} with $0\le \varphi_k < \pi$. Near $\theta = 0^+$ , the whole sum $f$ is equal to $f(\theta)=\sum_k |c_k|\sin(\theta +\varphi_k)$ which we write as $R\sin(\theta + \psi)$.

Each function $\sin(\theta + \varphi_k)$ changes sign twice in $(0, 2\pi]$, once when $\theta = t_k = \pi-\varphi_k$ and once when $\theta=t_k'=2\pi-\varphi_k$.

Moving $\theta$ from left to right, when we cross the value $t_k$, we add $-2|c_k|\sin(\theta + \varphi_k)$ to $f$, which leads to \begin{equation} R\sin(\theta+\psi) - 2 |c_k|\sin(\theta + \varphi_k) \quad \longrightarrow \quad R'\sin(\theta+\psi') \end{equation} Similarly, when we cross the value $t_k'$, we add $2|c_k|\sin(\theta + \varphi_k)$ to $f$.

Each addition to $f$ is a $O(1)$ step and the computation of the maximum of $f=R\sin(\theta+\psi)$ on the current interval is also a $O(1)$ step, thus the whole algorithm has $O(n)$ complexity, or $O(n \log n)$ if you count the initial step to sort the numbers $t_k$.