Finding maximum

First, let's reference the "brute force" force algorithm. This will have some issues, that I will call out below, but it is a correct solution.

struct Result
{
    size_t i;
    size_t j;
    int64_t value;
};

Result findBestBruteForce(const vector<int>& a)
{
    size_t besti = 0;
    size_t bestj = 0;
    int64_t bestvalue = INT64_MIN;

    for (size_t i = 0; i < a.size(); i++)
    {
        for (size_t j = i + 1; j < a.size(); j++)
        {
            // do the math in 64-bit space to avoid overflow
            int64_t value = (a[i] + (int64_t)a[j]) * (j - i);
            if (value > bestvalue)
            {
                bestvalue = value;
                besti = i;
                bestj = j;
            }
        }
    }
    return { besti, bestj, bestvalue };
}

The problem with the above code is that it runs at O(N²). Or more precisely, for the the N iterations of the outer for-loop (where i goes from 0 to N), there are an average of N/2 iterations on the inner for-loop. If N is small, this isn't a problem.

On my PC, with full optimizations turned on. When is N under 20000, the run time is less than a second. Once N approaches 100000, it takes several seconds to process the 5 billion iterations. Let's just go with a "billion operations per second" as an expected rate. If N were to 1000000, the maximum as the OP outlined, it would probably take 500 seconds. Such is the nature of a N-squared algorithm.

So how can we speed it up? Here's an interesting observation. Let's say our array was this:

10 5 4 15 13 100 101 6

On the first iteration of the outer loop above, where i=0, we'd be computing this on each iteration of the inner loop:

for each j: (a[0]+a[j])(j-0)
for each j: (10+a[j])(j-0)
for each j: [15*1, 14*2, 25*3, 23*4, 1000*5, 1010*6, 16*6]
         =  [15,     28,   75,   92,   5000,   6060,   96]

Hence, for when i=0, a[i] = 15 and the largest value computed from that set is 6060.

Since A[0] is 15, and we're tracking a current "best" value, there's no incentive to iterate all the values again for i=1 since a[1]==14 is less than 15. There's no j index that would compute a value of (a[1]+a[j])*(j-1) larger than what's already been found. Because (14+a[j])*(j-1) will always be less than (15+a[j])*(j-1). (Assumes all values in the array are non-negative).

So to generalize, the outer loop can skip over any index of i where A[best_i] > A[i]. And that's a real simple alteration to our above code:

Result findBestOptimized(const std::vector<int>& a)
{

    if (a.size() < 2)
    {
        return {0,0,INT64_MIN};
    }

    size_t besti = 0;
    size_t bestj = 0;
    int64_t bestvalue = INT64_MIN;
    int minimum = INT_MIN;

    for (size_t i = 0; i < a.size(); i++)
    {
        if (a[i] <= minimum)
        {
            continue;
        }

        for (size_t j = i + 1; j < a.size(); j++)
        {
            int64_t value = (a[i] + (int64_t)a[j]) * (j - i);
            if (value > bestvalue)
            {
                bestvalue = value;
                besti = i;
                bestj = j;
                minimum = a[i];
            }
        }
    }

    return { besti, bestj, bestvalue };
}

Above, we introduce a minimum value for A[i] to be before considering doing the full inner loop enumeration.

I benchmarked this with build optimizations on. On a random array of a million items, it runs in under a second.

But wait... there's another optimization!

If the inner loop fails to find an index j such that value > bestvalue, then we already know that the current A[i] is greater than minimum. Hence, we can increment minimum to A[i] regardless at the end of the inner loop.

Now, I'll present the final solution:

Result findBestOptimizedEvenMore(const std::vector<int>& a)
{

    if (a.size() < 2)
    {
        return { 0,0,INT64_MIN };
    }

    size_t besti = 0;
    size_t bestj = 0;
    int64_t bestvalue = INT64_MIN;
    int minimum = INT_MIN;

    for (size_t i = 0; i < a.size(); i++)
    {
        if (a[i] <= minimum)
        {
            continue;
        }

        for (size_t j = i + 1; j < a.size(); j++)
        {
            int64_t value = (a[i] + (int64_t)a[j]) * (j - i);
            if (value > bestvalue)
            {
                bestvalue = value;
                besti = i;
                bestj = j;
            }
        }
        minimum = a[i]; // since we know a[i] > minimum, we can do this

    }

    return { besti, bestj, bestvalue };
}

I benchmarked the above solution on different array sizes from N=100 to N=1000000. It does all iterations in under 25 milliseconds.

In the above solution, there's likely a worst case runtime of O(N²) again when all the items in the array are in ascending order. But I believe the average case should be on the order of O(N lg N) or better. I'll do some more analysis later if anyone is interested.

Note: Some notation for variables and the Result class in the code have been copied from @selbie's excellent answer.

Here's another O(n^2) worst-case solution with (likely provable) O(n) expected performance on random permutations and room for optimization.

Suppose [i, j] are our array bounds for an optimal pair. By the problem definition, this means all elements left of i must be strictly less than A[i], and all elements right of j must be strictly less than A[j].

This means we can compute the left-maxima of A: all elements strictly greater than all previous elements, as well as the right-maxima of A. Then, we only need to consider left endpoints from the left-maxima and right endpoints from the right-maxima.

I don't know the expectation of the product of the sizes of left and right maxima sets, but we can get an upper bound. The size of left maxima is at most the size of the longest increasing subsequence (LIS) of A. The right maxima are at most the size of the longest decreasing subsequence. These aren't independent, but I'm taking as an (unproven) assumption that the LIS and LDS lengths are inversely correlated with each other for random permutations. The right-maxima must start after the left-maxima end, so this seems like a safe assumption.

The length of the LIS for random permutations follows the Tracy-Widom distribution, so it has mean sqrt(2N) and standard deviation N^(-1/6). The expected square of the size is therefore 2N + 1/(N^1/3) so ~2N. This isn't exactly the proof we wanted, since you'd need to sum over the partial density function to be rigorous, but the LIS is already an upper bound on the left-maxima size, so I think the conclusion is still true.

C++ code (Result class and some variable names taken from selbie's post, as mentioned):

struct Result
{
    size_t i;
    size_t j;
    int64_t value;
};


Result find_best_sum_size_product(const std::vector<int>& nums)
{   
    /* Given: list of positive integers nums
       Returns: Tuple with (best_i, best_j, best_product)
       where best_i and best_j maximize the product
      (nums[i]+nums[j])*(j-i) over 0 <= i < j < n

       Runtime: O(n^2) worst case,
       O(n) average on random permutations.
    */
    
    int n = nums.size();
    
    if (n < 2)
    {
        return {0,0,INT64_MIN};
    }
    
    
    std::vector<int> left_maxima_indices;
    left_maxima_indices.push_back(0);
    
    for (int i = 1; i < n; i++){
        if (nums.at(i) > nums.at(left_maxima_indices.back())) {
            left_maxima_indices.push_back(i);
        }
    }
    
    std::vector<int> right_maxima_indices;
    right_maxima_indices.push_back(n-1);
    
    for (int i = n-1; i >= 0; i--){
        if (nums.at(i) > nums.at(right_maxima_indices.back())) {
            right_maxima_indices.push_back(i);
        }
    }
    
    size_t best_i = 0;
    size_t best_j = 0;
    int64_t best_product = INT64_MIN;
    int i = 0;
    int j = 0;

    for (size_t left_idx = 0;
         left_idx < left_maxima_indices.size();
         left_idx++)
    {
        i = left_maxima_indices.at(left_idx);

        for (size_t right_idx = 0;
            right_idx < right_maxima_indices.size();
            right_idx++)
        {   
            j = right_maxima_indices.at(right_idx);
            if (i == j) continue;
            
            int64_t value = (nums.at(i) + (int64_t)nums.at(j)) * (j - i);
            if (value > best_product)
            {
                best_product = value;
                best_i = i;
                best_j = j;
            }
        }
    }

    return { best_i, best_j, best_product };
}

I started from the two excellent answers by @selbie and @kcsquared.

Their solutions gave impressive results for random inputs. What was not clear is the worst case behavior.

What sequence would correspsond to the worst case?

I finally found a critial sequence for these two answers, a triangle sequence: this sequence slightly increases up to a max, and then slightly decrease. With such a sequence and n=10^5 for example, these answers take more than 10s.

My solutions starts from @selbie solution and add two improvements:

• I add @kcsquared's trick: on the right (of j), they can be only lower elements

• When considering a new left element a[i], it is useless to start from i + 1 to get the second element. We can start from the current best_j

With these tricks, I was able to improve the two posted answer performances a little bit. However, it still fails to solve the triangle sequence issue: about 10s for n = 10^5.

#include <iostream>
#include <vector>
#include <string>
#include <cstdlib>
#include <ctime>
#include <chrono>

struct Result {
    size_t i;
    size_t j;
    int64_t value;
};

void print (const Result& res, const std::string& prefix = "") {
    std::cout << prefix;
    std::cout << "(" << res.i << ", " << res.j << ") -> " << res.value << std::endl;
}

Result findBest(const std::vector<int>& a) {
    if (a.size() < 2) {
        return { 0, 0, INT64_MIN };
    }
    int n = a.size();  
    std::vector<int> next_max(n, -1);
    int current_max = n-1;
    for (int i = n-1; i >= 0; --i) {
        if (a[i] > a[current_max]) {
            current_max = i;
        }
        next_max[i] = current_max;
    }
    
    size_t besti = 0;
    size_t bestj = 0;
    int64_t bestvalue = INT64_MIN;
    int minimum = INT_MIN;
    
    for (size_t i = 0; i < a.size(); i++) {
        if (a[i] <= minimum) {
            continue;
        }
        minimum = a[i];
        size_t jmin = (bestj > i) ? bestj : i+1;
        for (size_t j = jmin; j < a.size(); j++) {
            j = next_max[j];
            value = (a[i] + (int64_t)a[j]) * (j - i);
            if (value > bestvalue) {
                bestvalue = value;
                besti = i;
                bestj = j;
            } 
        }
    }
    return { besti, bestj, bestvalue };
}

int main() {
    int n = 1000000;
    int vmax = 100000000;
       
    std::vector<int> A (n);
    std::srand(std::time(0));
    for (int i = 0; i < n; ++i) {
        A[i] = rand() % vmax + 1;
    }
    
    std::cout << "n = " << n << std::endl;
    auto t0 = std::chrono::high_resolution_clock::now();
    auto res = findBest (A);
    auto t1 = std::chrono::high_resolution_clock::now();
    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0).count();
    print (res, "Random: ");
    std::cout << "time = " << duration/1000 << " ms" << std::endl;
    
    int i_max = n/2;
    for (int i = 0; i < i_max; ++i) A[i] = i+1;
    A[i_max] = 10 * i_max;
    for (int i = i_max+1; i < n; ++i) {
        A[i] = 2*i_max - i;
    }
    t0 = std::chrono::high_resolution_clock::now();
    res = findBest (A);
    t1 = std::chrono::high_resolution_clock::now();
    duration = std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0).count();
    print (res, "Triangle sequence: ");
    std::cout << "time = " << duration/1000 << " ms" << std::endl;
    
    return 0;
}