Finding maximum

Question

Given an integer n and array a. Finding maximum of (a[i]+a[j])*(j-i) with 1<=i<=n-1 and i+1<=j<=n

Example:

Input

5
1 3 2 5 4

Output

21

Explanation :With i=2 and j=5, we have the maximum of (a[i]+a[j])*(j-i) is (3+4)*(5-2)=21

Constraints:

n<=10^6
a[i]>0 with 1<=i<=n

I can solve this problem with n<=10^4, but what should I do if n is too large, like the constraints?

Answer 1

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 };
}

Answer 2

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;
}