Optimizar funciones unimodales

Introducción

En programación competitiva, a menudo se nos pide hallar el valor mínimo o máximo que podemos alcanzar bajo ciertas condiciones. Más específicamente, a menudo se nos pide formular alguna función $f(x)$ y hallar el valor máximo o mínimo de $f(x)$.

Si nuestra función resulta ser unimodal, podemos usar o bien búsqueda ternaria o bien búsqueda binaria para hallar la moda de la función de forma eficiente.

Definiciones

Una función es unimodal si se cumple una de las siguientes condiciones:

1. La función crece de forma estricta primero, luego alcanza un máximo, luego decrece de forma estricta.
2. La función decrece de forma estricta primero, luego alcanza un mínimo, luego crece de forma estricta.

Asumamos que queremos hallar el mínimo local de nuestra función. Entonces, una función es convexa si las diferencias consecutivas son no decrecientes. En términos de cálculo, si la derivada de la función es no decreciente, la función es convexa.

Nótese que esto implica que una función convexa también es una función unimodal.

Búsqueda ternaria

Por ahora, asumamos que queremos hallar el mínimo de nuestra función, y que nuestra función es unimodal. Esto significa que se cumple la segunda condición: nuestra función decrece de forma estricta primero, alcanza su mínimo y luego crece de forma estricta.

La idea detrás de la búsqueda ternaria es un tanto similar a la búsqueda binaria. En cada paso de nuestro algoritmo, queremos recortar una porción grande de nuestro espacio de búsqueda para dejar de considerarla.

Consideremos dos puntos $m_1$ y $m_2$, donde $l < m_1 < m_2 < r$. Aquí, $l$ y $r$ son los extremos de nuestro espacio de búsqueda actual. Esto termina dividiendo nuestro espacio de búsqueda en tres secciones.

Ahora, hacemos un poco de análisis por casos para ver si podemos eliminar alguna sección.

1. $f(m_1) > f(m_2)$ significa que podemos eliminar $[l, m_1)$
2. $f(m_1) < f(m_2)$ significa que podemos eliminar $(m_2, r]$
3. $f(m_1) = f(m_2)$ significa que podemos eliminar $[l, m_1)$ y $(m_2, r]$

Típicamente, tratamos el caso 3 como análogo a los casos 1 y 2.

Como mini-demostración de por qué la forma en que manejamos el caso 1 es correcta, consideremos lo siguiente:

Si $f(m_1) > f(m_2)$, entonces sabemos estos dos hechos:

• $f(m_1)$ no es el mínimo
• $f(m_1)$ no puede estar en el lado estrictamente creciente de la función Respecto al segundo punto, si $f(m_1) > f(m_2)$ y $m_1$ estuviera en el lado estrictamente creciente de la función, eso implicaría $m_1 > m_2$, lo cual es una contradicción.

Así, podemos concluir que $m_1$ permanece en la sección estrictamente decreciente de nuestra función. Por esto, podemos concluir que nuestro mínimo no estará en el rango $[l, m_1)$, así que podemos eliminar esta sección de nuestro espacio de búsqueda.

Ahora que sabemos cómo eliminar secciones de nuestro espacio de búsqueda de forma estratégica, queda la pregunta de cómo elegir los mejores valores de $m_1$ y $m_2$. La forma óptima es dividir nuestro espacio de búsqueda en tercios, así que tenemos:

Implementación

La complejidad temporal de la búsqueda ternaria sigue la siguiente recurrencia:

Por el teorema maestro , esto es $\mathcal{O}(\log{n})$ si trabajamos con enteros.

Si hacemos búsqueda ternaria con un margen de error fijo $\epsilon$, la complejidad temporal pasa a ser $\mathcal{O}\left( \log\left( \frac{n}{\epsilon} \right) \right)$.

template <typename F> double find_min(double l, double r, double eps, const F &f) {
	while (r - l > eps) {
		double m1 = l + (r - l) / 3;
		double m2 = r - (r - l) / 3;
		f(m1) > f(m2) ? l = m1 : r = m2;
	}
 
	return l;
}

def find_min(l: float, r: float, epsilon: float, f) -> float:
	while r - l > epsilon:
		m1 = l + (r - l) / 3
		m2 = r - (r - l) / 3
		if f(m1) > f(m2):
			l = m1
		else:
			r = m2
	return l

Si nuestra función solo toma enteros, entonces la implementación se ve un poco distinta.

template <typename F> int find_min(int l, int r, const F &f) {
	while (r - l > 3) {
		int m1 = l + (r - l) / 3;
		int m2 = r - (r - l) / 3;
		f(m1) > f(m2) ? l = m1 : r = m2;
	}
 
	int res = l;
	for (int i = l + 1; i <= r; i++) {
		if (f(i) < f(res)) { res = i; }
	}
 
	return res;
}

def find_min(l: int, r: int, f) -> int:
	while l < r:
		m1 = l + (r - l) // 3
		m2 = r - (r - l) // 3
		if f(m1) > f(m2):
			l = m1 + 1
		else:
			r = m2 - 1
	return l

Búsqueda binaria

De hecho, cuando nuestra función solo toma enteros, la búsqueda binaria suele ser una mejor opción, porque es más corta y requiere menos llamadas a nuestra función $f(x)$.

Sin embargo, si trabajamos con punto flotante, la búsqueda binaria puede ser un poco más problemática por la pérdida de precisión. Por esta razón, la búsqueda ternaria a veces sigue siendo preferible.

template <typename F> int find_min(int l, int r, const F &f) {
	while (l < r) {
		int m = (l + r) / 2;
		f(m) < f(m + 1) ? r = m : l = m + 1;
	}
 
	return l;
}

def find_min(l: int, r: int, f) -> int:
	while l < r:
		m = (l + r) // 2
		if f(m) > f(m + 1):
			l = m + 1
		else:
			r = m
	return l

La forma en que funciona el enfoque de búsqueda binaria es que halla el primer punto donde la función se vuelve estrictamente creciente. Esto funciona porque asumimos que nuestra función toma la forma de decrecer de forma estricta, alcanzar su mínimo y luego crecer de forma estricta. Así, el primer punto donde nuestra función se vuelve estrictamente creciente es el punto de inflexión de nuestra función, que es el mínimo.

Ejemplo - Haybale Distribution

Explicación

Definamos una función $f(y)$ que evalúa el costo de transportar todos los fardos de heno si los entregamos en el punto $y$. ¿Podemos mostrar que esta función es unimodal o convexa?

Un hecho importante al mostrar que una función es convexa es que la suma de funciones convexas es convexa.

Aquí, $f(y)$ es la suma de todas las funciones de costo individuales de cada granero. La función de costo de un granero dado es convexa, ¡así que resulta que $f(y)$ también es convexa!

Dado que una función convexa es unimodal, esto nos permite hacer búsqueda ternaria sobre ella. Así, podemos hallar el mínimo en $\mathcal{O}(\log{N})$ por consulta.

Implementación

Complejidad temporal: $\mathcal{O}((N + Q) \log{N})$

Nótese que la implementación de abajo usa búsqueda binaria en su lugar, porque es un poco más fácil de implementar.

#include <bits/stdc++.h>
 
using ll = long long;
 
int main() {
	int n;
	std::cin >> n;
	std::vector<int> x(n);
	for (int &i : x) { std::cin >> i; }
 
	std::sort(x.begin(), x.end());
 
	// compute prefix sums on x[i] values
	std::vector<ll> pref(n + 1);
	for (int i = 1; i <= n; i++) { pref[i] = pref[i - 1] + x[i - 1]; }
 
	/** @return sum of x[i] on range [l, r] */
	const auto pref_sum = [&](int l, int r) -> ll { return pref[r + 1] - pref[l]; };
 
	// only keep track of distinct x[i]
	std::vector<int> unique_x = {0};
	for (int i = 1; i < n; i++) {
		if (x[i] != x[i - 1]) { unique_x.push_back(i); }
	}
 
	int query_num;
	std::cin >> query_num;
	for (int q = 0; q < query_num; q++) {
		int a, b;
		std::cin >> a >> b;
 
		/** @return cost of delivering hay at x[unique_x[idx]] */
		const auto cost = [&](int idx) -> ll {
			int loc = unique_x[idx];
			ll left_dists = 1ll * x[loc] * loc - pref_sum(0, loc - 1);
			ll right_dists = pref_sum(loc, n - 1) - 1ll * x[loc] * (n - loc);
			return left_dists * a + right_dists * b;
		};
 
		/*
		 * We use binary search to find the minimum (instead of ternary).
		 * In this case, we find the first point where the function becomes increasing.
		 *
		 * To optimize each query to O(log n), we choose to only consider
		 * unique x[i] values to deliver our haybales at.
		 */
		int l = 0, r = unique_x.size() - 1;
		while (l < r) {
			int m = (l + r) / 2;
			cost(m) < cost(m + 1) ? r = m : l = m + 1;
		}
 
		std::cout << cost(l) << '\n';
	}
}

Problemas