SIMARRAY - Editorial

PROBLEM LINK:

Practice
Div-2 Contest
Div-1 Contest

Author: Sundar
Tester: Riley Borgard
Editorialist: Jatin Yadav

DIFFICULTY:

MEDIUM

PREREQUISITES:

Quadratic Functions

PROBLEM:

You are given two arrays of positive integers, A and B, both of length N. Find the minimum possible value of \displaystyle \sum_{i=1}^{N} (A_i - x_i B_i)^2, over all possible N-tuples of real numbers (x_1, x_2, \ldots x_N) satisfying x_1 \ge x_2 \ge \cdots \ge x_N

EXPLANATION:

For each 1 \le L \le R \le N, define F_{[L,R]}(x) = \sum_{i=L}^{R} (A_i - B_ix)^2 = Q_{[L,R]} x^2 - 2 P_{[L,R]} x + C_{[L,R]}, where Q_{[L,R]} is the sum of B_i^2, P_{[L, R]} is the sum of A_iB_i and C_{[L,R]} is the sum of A_i^2 in the range [L,R] respectively.

It is minimized at t_{[L,R]} = \dfrac{P_{[L,R]}}{Q_{[L,R]}}. For all x < t_{[L,R]}, it is decreasing and for all x > t_{[L,R]}, it is increasing. Let the minimum value be V_{[L,R]}

Define x_0 = \infty, x_{N+1} = -\infty. Lets prove the following statement:

Lemma 1 : Lets partition the optimal solution into maximal segments of equal values. Consider any such segment [L,R]. Then, x_i = t_{[L,R]} for all i in [L,R].

Proof : Assume that all x_i for i in range [L,R] were equal to u \ne t_{[L,R]}. Let r = \min(|u - t_{[L,R]}|, \min(x_{L-1} - u, u - x_{R+1})). Clearly, r > 0. If u < t_{[L,R]} we can change all values to u+r, else we can change all values to u-r to improve the objective, while maintaining the non increasing property.

So, any optimal solution can be thought of as segments [1=L_1,R_1], [L_2=R_1+1,R_2] \ldots [L_{k} = R_{k-1}+1,R_k = N] with the segment [L_i,R_i] having a value x_i = t_{[L_i,R_i]}, such that x_1\geq x_2 \ldots \geq x_k.


Subtasks 1 and 2

There are exactly 2^{N-1} ways of partitioning the array into segments (for each index i > 1, choose whether a new segment starts at i). So, we can brute-force over these partitions and find the minimum objective over all the valid partitions. This passes subtasks 1 and 2.


Subtask 3

For subtask 3, let’s define dp[L][R] to be the minimum value of the sums of deviation over first R indices, if the last segment was [L,R]. Then, let’s iterate over the left end L' of the second last segment, such that t_{[L',L-1]} \geq t_{[L,R]}. Over all such L', we have to minimize dp[L'][L-1] + V_{[L,R]}. The total time complexity would be O(N^3).


Subtask 4

We need to minimize dp[L'][L-1] over all L' with t_{[L',L-1]} \geq t_{[L,R]}. Let’s sort the intervals [L',L-1] according to t_{[L,L-1]}. Then, we need the minimum of a suffix. The suffix minimums can be precomputed. We don’t have to sort these segments ending at L-1 for every R. We can do this once we have processed all segments with right end L-1. The time complexity would be O(N^2 \log(N))


Subtask 5

Let’s consider a more general problem. Say we have m quadratic functions H_1(x), H_2(x), \ldots H_m(x) each having a positive coefficient of x^2. Let the respective minimizing points be v_1, v_2, \ldots, v_m. We want to minimize \sum H_i (x_i) over all tuples x with x_1 \geq x_2 \ldots \geq x_m .

Lemma 2: Consider some i with v_i < v_{i+1}. Then, in the optimal solution, x_i = x_{i+1}

Proof : Assume that in the optimal solution x_i > x_{i+1}. If x_i > v_i, we can decrement x_i by a small enough value \epsilon improving the objective without violating the non increasing constraint. Similarly, if x_{i+1} < v_{i+1}, we can increment x_{i+1} by a small enough \epsilon to improve the objective. Therefore, x_i \le v_i and x_{i+1} \geq v_{i+1}. But then, v_{i+1} \le x_{i+1} < x_i \leq v_i, a contradiction.

Using this, we can get a simple O(N^2) algorithm. If there exists an index with v_i < v_{i+1}, remove H_i(x), H_{i+1}(x) from the sequence and add H_i(x) + H_{i+1}(x).

Once v_i \geq v_{i+1} for all i, the solution x_i = v_i for all i is valid and clearly optimal.


Subtask 6

To improve the above algorithm, we can just keep a stack of quadratic functions. For each i, we first add H_i(x) to the stack and then keep on merging the top two elements of the stack while their minimizing points are out of order. This leads to a simple O(N) algorithm.

SOLUTIONS:

Setter's Solution
#include <bits/stdc++.h>
using namespace std;
 
#define ll long long
#define pii pair<int, int>
#define all(c) ((c).begin()), ((c).end())
#define sz(x) ((int)(x).size())
 
#ifdef LOCAL
#include <print.h>
#else
#define trace(...)
#endif
 
long long readInt(long long l,long long r,char endd){
	long long x=0;
	int cnt=0;
	int fi=-1;
	bool is_neg=false;
	while(true){
		char g=getchar();
		if(g=='-'){
			assert(fi==-1);
			is_neg=true;
			continue;
		}
		if('0'<=g && g<='9'){
			x*=10;
			x+=g-'0';
			if(cnt==0){
				fi=g-'0';
			}
			cnt++;
			assert(fi!=0 || cnt==1);
			assert(fi!=0 || is_neg==false);
 
			assert(!(cnt>19 || ( cnt==19 && fi>1) ));
		} else if(g==endd){
			if(is_neg){
				x= -x;
			}
			if(!(l<=x && x<=r))cerr<<l<<"<="<<x<<"<="<<r<<endl;
			assert(l<=x && x<=r);
			return x;
		} else {
			assert(false);
		}
	}
}
string readString(int l,int r,char endd, char minc = 'a', char maxc = 'z'){
	string ret="";
	int cnt=0;
	while(true){
		char g=getchar();
		assert(g!=-1);
		if(g==endd){
			break;
		}
		assert(g >= minc && g <= maxc);
		cnt++;
		ret+=g;
	}
	assert(l<=cnt && cnt<=r);
	return ret;
}
long long readIntSp(long long l,long long r){
	return readInt(l,r,' ');
}
long long readIntLn(long long l,long long r){
	return readInt(l,r,'\n');
}
string readStringLn(int l,int r, char minc = 'a', char maxc = 'z'){
	return readString(l,r,'\n', minc, maxc);
}
string readStringSp(int l,int r, char minc = 'a', char maxc = 'z'){
	return readString(l,r,' ', minc, maxc);
}
template<class T>
vector<T> readVector(int n, long long l, long long r){
    vector<T> ret(n);
    for(int i = 0; i < n; i++){
        ret[i] = i == n - 1 ? readIntLn(l, r) : readIntSp(l, r);
    }
    return ret;
}
 
struct Rational
{
    ll num, den;
    int cnt;
    Rational(ll n, ll d, int cnt) : num(n), den(d), cnt(cnt){}
    bool operator<(const Rational &r) const{
        return num * __int128_t(r.den) < den * __int128_t(r.num);
    };
    Rational operator + (const Rational & r) const{
        return Rational(num + r.num, den + r.den, cnt + r.cnt);
    }
};
 
const int MAXN = 500000;
int main(){
    int t = readIntLn(1, MAXN);
    int sn = 0;
    while(t--){
        int n = readIntLn(2, MAXN); 
        sn += n;
        assert(sn <= MAXN);
        vector<ll> a = readVector<ll>(n, 1, 1000), b = readVector<ll>(n, 1, 1000);
        double ans = 0;
 
        for(ll & x : a){
            ans += x * x;
        }
 
        vector<Rational> stk;
 
        for(int i = 0; i < n; i++){
            a[i] *= b[i];
            b[i] *= b[i];
            // b * x^2 - 2 * a * x
            stk.push_back(Rational(a[i], b[i], 1));
            while(stk.size() >= 2){
                Rational x = stk.back(); stk.pop_back();
                Rational y = stk.back(); stk.pop_back();
                if(y < x){
                    stk.push_back(x + y);
                } else{
                    stk.push_back(y);
                    stk.push_back(x);
                    break;
                }
            }
        }
        int j = 0;
        for(int i = 0; i < stk.size(); i++){
            double x = stk[i].num / (double) stk[i].den;
            for(int k = 0; k < stk[i].cnt; k++, j++){
                ans += x * (b[j] * x - 2 * a[j]);
            }
        }
        trace(stk.size());
        cout << setprecision(20) << fixed << ans << endl;
    }
}
Tester's Solution
#include <bits/stdc++.h>
 
#define ll long long
#define sz(x) ((int) (x).size())
#define all(x) (x).begin(), (x).end()
#define vi vector<int>
#define pii pair<int, int>
#define rep(i, a, b) for(int i = (a); i < (b); i++)
using namespace std;
template<typename T>
using minpq = priority_queue<T, vector<T>, greater<T>>;
 
using ftype = long double;
 
struct node {
    ftype ab, bb;
    int n;
    node(ftype a, ftype b) {
        n = 1;
        ab = a * b;
        bb = b * b;
    }
    node(ftype ab, ftype bb, int n) : ab(ab), bb(bb), n(n) {}
    node operator+(const node &o) const {
        return node(ab + o.ab, bb + o.bb, n + o.n);
    }
    ftype getx() const {
        return ab / bb;
    }
};
 
void solve() {
    int n;
    cin >> n;
    vector<ftype> a(n), b(n);
    rep(i, 0, n) cin >> a[i];
    rep(i, 0, n) cin >> b[i];
    vector<node> st;
    rep(i, 0, n) {
        node x(a[i], b[i]);
        while(!st.empty() && st.back().getx() < x.getx()) {
            x = x + st.back();
            st.pop_back();
        }
        st.push_back(x);
    }
    vector<ftype> r;
    for(node &x : st) {
        ftype val = x.getx();
        rep(i, 0, x.n) r.push_back(val);
    }
    ftype ans = 0;
    rep(i, 0, n) ans += (a[i] - r[i] * b[i]) * (a[i] - r[i] * b[i]);
    cout << ans << '\n';
}
 
int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);
    cout << fixed << setprecision(10);
    int te;
    cin >> te;
    while(te--) solve();
}