I think it is possible to solve this problem in $\mathcal O(n \log n)$. Note that this is an improvement to this answer, which previously did not handle that $i \ne j$. See below for how this is added. Also note that what described below is a version of the Convex Hull Trick.
First, observe that the $\max\{a_i b_j, a_j b_i\}$ is unnecessary: It can be replaced by $a_i b_j$, since we get the other result by swapping $i, j$.
This leaves us with finding this:
$$\max_{i, j} \quad a_ib_i + a_jb_j + a_ib_j $$
Our approach will now be to iterate over $i$ and find the largest $j$ in $\mathcal O(\log n)$ time. Note that for this, it suffices to, given $c$, find $j$ maximizing this:
$$\max_{j} \quad (a_j + c) b_j $$
(Set $c := a_i$ and observe that $a_i b_i$ is not dependent on $j$)
Now notice that this defines a collection of functions $f_j(x) = a_j b_j + x b_j$, which are all linear (plus offset). We now want to find which of these functions is the largest, given $x$. We now build a data structure that can solve this quickly (i.e. in $\log n$).
First notice that for functions with the same $b_j$, the one with the highest $a_j$ is always higher. Geometrically, these are parallel lines, but those with higher $a_j$ are higher. We thus ignore all such dominated functions.
Then notice that a function with higher $b_j$ eventually dominates those with smaller $b_j$. We will now build a ordered list, where the function that dominates for $x \to \infty$ comes first, the ones that dominates to the left of it comes second, etc. The list will only contain those functions that dominate somewhere. Graphically, the list would contain [black, blue, red, green] , since the functions in this image below would dominate in that order from right to left. orange does not appear, since it dominates nowhere:
So we start with the function with the highest $b_j$ and add it to the list. The list is ordered such that the functions dominating more left come last, and we add from the back¹. Not that each function dominates in a unique interval. We also store the right end of this interval (which is $+\infty$ for the first function).
We now add the other functions, in order of decreasing $b_j$. (This means we sort the arrays once, which is $\mathcal O(n \log n)$).
To add such a function $f$, compare it with the function $g$ currently dominating the most on the left (which is the last function in the list). There are several options, depending on where $g$ starts dominating $f$ (which we can compute easily with some linear algebra):
- $f$ dominates $g$ on some interval, but not on the entire interval where $g$ dominates. It still leaves space for $g$ to dominate somewhere. In this case, we add $f$ to the back of the list.
- $f$ dominates $g$ everywhere. In this case, we remove $g$ (we drop $g$), and continue with the rest of the list.
- $f$ dominates $g$ nowhere: Technically, the point of domination is $<0$. As an optimization, we can choose to do leave the list unchanged (but it does not really matter, the algorithm should still work for numbers $<0$).
- As a special case, we can also do the pruning of smaller functions with equal $b_j$ here: If two functions have same $b_j$, process those with higher $a_j$ first, and drop the others here since they never dominate.
The runtime analysis here is easy, once you see that each step either adds a function directly, or drops a function, where each function is only dropped once. Thus by amortized analysis, where each function in the list has a potential of $2$ (and adding thus requires a potential of $3$), we see that this is works in $\mathcal O(n)$.
We now have a sorted list that tells us where each function dominates. Thus, we can find the highest function for given $x$ in $\mathcal O(\log n)$ by binary search.
This means we're done. We can find for each $i$ the highest $j$ by doing the above binary search with $x := c := a_i$. We have $n$ iterations of the outer loop, each taking $\mathcal O(\log n)$. We also have $\mathcal O(n \log n) + \mathcal O(n)$ (for sorting and building the function sorted list) of setup time. Together, that's $\mathcal O(n \log n)$.
To handle $i \ne j$, we must additionally keep track of the second-highest function each. This seems to actually not be that hard, the trick is to maintain a similar data structure for the second-highest function.
When adding a function $f$, all the functions that we dropped are now the second-highest functions on the interval left of the intersection point of $f$ with $g$ (the function it non-trivially intersects with). But first, we notice that right of the intersection point, $f$ is now a candidate for the second-highest function, that is also higher than all second-highest functions left of the intersection point. So we can add it to the second-highest-function lookup structure (using a regular insert), and then force-add all the things we dropped, that were to the left of the intersection point. What is tricky here is that there are various edge-cases, and that the data structure is no longer as well-behaved (since when going right, functions can become less steep).
This idea actually generalized to also keeping track of the $m$-th highest list. The runtime complexity for an insert should then be $\mathcal O(n m)$ using amortized analysis. For us, $m=2$ is sufficient, making it linear.
I've implemented this, and tested it a bit, and it seems to work. But I am sure there are bugs somewhere, and it's quite possible that this approach is wrong in general. But have a look yourself here: https://github.com/JoJoDeveloping/some-optimization-problem
Unfortunately, I don't see how to extend this to handling $i \neq j$. Perhaps you can somehow save the second-hightest function for each point? But I don't know how to build this lookup structure efficiently. The second-highest function is also way less well-behaved (observe that blue is the second-highest function for two disjoint intervals in the picture above).
1: This allows us to use an array. Alternatively, we can use a linked list and convert it to an array later.
