LP

정의

제약조건과 목적함수가 모두 "선형"인 계의 최적화방법에 대한 이론이다.

요리조리 식정리를 하면 반드시 다음과 같은 형태로 만들 수 있다고 한다.

$Ax \leq b \\ x \geq 0 \\ \text{maximize} \ c^Tx$

식정리 cheat sheet

  • min ojb → -(max -obj)
  • a1x1+a2x2 ≥ b → a1x1+a2x2-x3=b
  • 부등호 뒤집기 → 양변에 -1곱
  • non-neg 조건이 없는 x → x=x1-x2, x1,x2>=0으로 치환
  • 목적함수의 변수에 절댓값 → 조건x=x1-x2,목적x=x1+x2 (최적화해법1 책 p.116참고)(slope trick이 일종의 LP인 이유)
  • 그외 → 2차계획법(목적함수2차,조건1차)도 된다는듯.
  • 조건개수 > 변수개수 → 듀얼시키는게 이득
  • min(max(~)) or max(min(~))형태 목적함수도 max나 min하나로 풀어줄 수 있음(최적화해법1 p.118 참고). 혹은 안쪽LP를 듀얼시켜 maxmax(minmin)으로 맞춰줄 수도 있음(13091번 참고)
  • 어떤 변수가 다른 변수를 참조하는 조건이어도 그냥 다 쓰고나서, 변수왼쪽 상수오른쪽으로 정리하면 된다. (13897번 모델링 참고)

활용

정수LP같은거 포함하면 끝도 없고, 정확히 LP에만 속하는건 최단거리,네트워크유량,일부 그리디 정도?

최적화문제(LP뿐만 아니라 다른것도)에는 Dual이라는 개념이 있다. 게임이론에 minimax과 maximin, 그래프이론에 최소정점덮개와 최대독립집합 등 알게모르게 여기저기 낑겨있는 개념이다. LP의 Dual의 최적해=원래문제의 최적해 라는 좋은 성질이 존재할 뿐만 아니라 해당 시스템 자체에 대한 여러가지 정보를 관찰할 수도 있다. 어떤 문제들은 Primal상태에서 푸는게 도저히 불가능해보이는데 Dual시키면 풀리기도 하고, Simplex법도 Dual시키면 조건개수가 줄어 더 빠르게 동작할수 있다거나. 듀얼시키는 방법도 간단한 공식이라 알아두면 여러모로 좋다.

알고리즘

Simplex

최악의 경우 지수시간이 걸리는 데이터가 존재하지만, 많은 현실적 문제들은 그러한 경우가 거의 없으며 평균적으로 굉장히 빠르게 동작한다. 특히 PS에서는 문제특징에따라 오래걸리는 데이터를 만드는게 불가능한 경우도 있고 그래서 가끔 플로우문제를 심플렉스로 비비적거리면 뚫을수도 있다고(...) (억지로 메모리제한을 줄여 막을수는 있다)

대략적인 작동방식은 LP가 만드는 Hyper Space의 Polytope형태(더 엄밀하게는 Simplex형태)의 외부꼭짓점에서 그 근방꼭짓점중에 목적함수가 최소화(최대화)되는 쪽으로 이동하는걸 반복하는게 끝이다(이것의 수학적표현을 이해하는게 어려울뿐,,,). 이 알고리즘이 올바르게 작동하려면 볼록성이 보장되어야 하며, 그렇기때문에 선형시스템인 LP에 사용될 수 있다.

정리한 식의 b가 전부 non-negative이면 slack variable m개 추가해서 정규형으로 만들어줄 수 있고, 아닐때는 2Phase법이나 BigM법으로 정규형으로 변환해줄 수 있다. BigM은 수치적 안정성이 낮아 상용프로그램에서는 안쓴다고 하는데 PS범위에서는 별 문제 없는듯(지금까지는 M=10000이하로 다 풀수 있었다). bigM법이 이해하는데나 구현이나 훨씬 쉽기 때문에 추천하고, 2phase법은 구현이 많이 긴건 아닌데 이해도 잘 안되서 굳이 필요없을듯 하다.

Interior Point

최악의 경우에도 다항시간내에 풀림이 보장되는 알고리즘. 대충 듣기로는 타원어쩌구 저쩌구해서 점점 최적해 근방으로 이동한다는데 당연히 상수가 심플렉스에 비해 매우 크기 때문에 상당히 큰 LP문제에서나 Simplex보다 빠르다고. PS하면서는 볼일없을것같다.

구현

#pragma once
#include "core/base.h"
#include "misc/util.h"

//입력: Ax<=b, obj
//출력: maximize obj*x
//numeric stability is sensitive by M
template<class T=f64,int M>
void dualize(Arr<Arr<T>> &a,Arr<T> &b,Arr<T>& obj){
	int m=sz(a), n=sz(a[0]);
	transpose(a),swap(b,obj);
	for(int i=0;i<n;i++){
		for(auto& j:a[i])j=-j;
		b[i]=-b[i];
	}
	for(auto& i:obj)i=-i;
}
template<class T=f64,int M>
T simplex(Arr<Arr<T>>& a,Arr<T>& b,Arr<T>& obj){
	int m=sz(a),n=sz(a[0]),s=0;
	if(m>n){dualize<T,M>(a,b,obj);return -simplex<T,M>(a,b,obj);}
	func(void,elim,int r1,int r2,int c){//elim r2
		if(r1==r2){T x=a[r1][c]; for(auto& i:a[r1])i/=x;}
		else{
			T x=a[r2][c]/a[r1][c]; if(-eps<x&&x<eps)return;
			for(int i=0;i<n+s+m+2;i++)
				a[r2][i]-=x*a[r1][i];
		}
	};

	//make all b>=0
	Arr<char> geq(m);
	for(int i=0;i<m;i++)
		if(b[i]<0){
			for(auto& j:a[i])j=-j;
			for(auto& r:a)r.emplb(0);
			a[i][-1]=-1,b[i]=-b[i],geq[i]=true,s++;
		}

	//n vars, s slacks(-1), m slacks(1), 1 z, 1 b_value
	Arr<int> p(m);//행의 기본변수
	obj.resize(n+s+m+2);
	for(int i=0;i<m;i++)
		a[i].resize(n+s+m+2),a[i][p[i]=n+s+i]=1,a[i][-1]=b[i],obj[p[i]]=geq[i]?-M:0;
	
	//z=f(x) == z-f(x)=0
	for(auto &i:obj)i=-i;
	obj[-2]=1;
	a.emplb(obj);
	for(int i=0;i<m;i++)
		elim(i,m,p[i]);
	
	//now shape of a = (m+1)*(n+s+m+2)
	while(true){
		int ev=0,lvi=-1;
		for(int i=0;i<n+s+m;i++)
			ev=a[-1][ev]>a[-1][i]?i:ev;
		if(a[-1][ev]>-eps)break;
		for(int i=0;i<m;i++)
			if(a[i][ev]>eps and (!~lvi or a[i][-1]/a[i][ev]<a[lvi][-1]/a[lvi][ev]))
				lvi=i;
		if(!~lvi or (-eps<a[lvi][ev]&&a[lvi][ev]<eps)) throw "unbounded";
		for(int i=0;i<m+1;i++)elim(lvi,i,ev);
		p[lvi]=ev;
	}
	//역추적 잘 모르겠다. 가우스소거 필요한듯?
	//if(?) throw "infeasible"
	return a[-1][-1];
}

문제접근 팁

변수의 개수보다 조건의 개수가 실행속도에 영향이 더 크므로 그렇다면 Dual시켜 푸는게 훨씬 빠르다.

일반적으로 시간은 변수200개 조건200*200개 듀얼에서 1초정도, 메모리는 200MB정도 든다. 심플렉스 막으려고 작정하고 메모리를 줄여버린 문제들은 시간이 넉넉해도 도저히 풀수가 없다... 그럴땐 포기하고 정해로 풀자.

문제 특성상 몇몇조건을 커팅해서 제거할수 있을때는 무조건 제거해주자. ex) 16833번은 평면위의 점이라는 특징때문에 커팅하면 매우 빨라진다.

관련문제