연립 일차방정식을 푸는 가우스 방법¶
$m$개의 미지수를 가진 $n$개의 선형 대수 방정식으로 이루어진 연립 방정식(SLAE)이 주어집니다. 이 연립 방정식을 풀어야 합니다. 즉, 해가 없는지, 정확히 하나의 해를 가지는지, 아니면 무수히 많은 해를 가지는지 판별해야 합니다. 그리고 해가 적어도 하나 존재하는 경우, 그중 아무것이나 하나를 구해야 합니다.
형식적으로 이 문제는 다음과 같이 정의됩니다. 다음 연립 방정식을 푸십시오.
여기서 계수 $a_{ij}$($i$는 1부터 $n$까지, $j$는 1부터 $m$까지)와 $b_i$($i$는 1부터 $n$까지)는 알려져 있으며, 변수 $x_i$($i$는 1부터 $m$까지)는 미지수입니다.
이 문제는 간단한 행렬 표현으로도 나타낼 수 있습니다.
여기서 $A$는 계수 $a_{ij}$로 이루어진 $n \times m$ 크기의 행렬이고, $b$는 크기가 $n$인 열 벡터입니다.
이 글에서 소개하는 방법은 임의의 수 p를 법으로 하는 방정식을 푸는 데에도 사용할 수 있다는 점을 짚어 둘 만합니다. 즉, 다음과 같습니다.
가우스¶
엄밀히 말하면, 아래에서 설명하는 방법은 "가우스-조르당(Gauss-Jordan)" 또는 가우스-조르당 소거법이라고 불러야 합니다. 이는 1887년에 조르당이 기술한 가우스 방법의 변형이기 때문입니다.
개요¶
이 알고리즘은 각 방정식에 변수가 하나만 남을 때까지 각 방정식에서 변수를 sequential elimination하는 것입니다. $n = m$이면, 이는 행렬 $A$를 단위 행렬로 변환한 뒤, 해가 유일하며 계수 $b_i$와 같은 이 자명한 경우의 방정식을 푸는 것으로 생각할 수 있습니다.
가우스 소거법은 두 가지 간단한 변환에 기반합니다.
- 두 방정식을 서로 바꿀 수 있습니다
- 임의의 방정식은 그 행(0이 아닌 계수를 곱한 것)과 다른 몇몇 행(임의의 계수를 곱한 것)의 선형 결합으로 대체할 수 있습니다.
첫 번째 단계에서 가우스-조르당 알고리즘은 첫 번째 행을 $a_{11}$로 나눕니다. 그다음, 알고리즘은 첫 번째 열의 계수가 모두 0이 되도록 첫 번째 행을 나머지 행에 더합니다. 이를 위해 i번째 행에는 첫 번째 행에 $- a_{i1}$을 곱한 것을 더해야 합니다. 이 연산은 벡터 $b$에도 똑같이 수행해야 한다는 점에 유의하십시오. 어떤 의미에서는, 벡터 $b$가 행렬 $A$의 $m+1$번째 열인 것처럼 동작합니다.
그 결과, 첫 번째 단계가 끝나면 행렬 $A$의 첫 번째 열은 첫 번째 행에서는 $1$, 다른 행에서는 $0$으로 이루어집니다.
마찬가지로 알고리즘의 두 번째 단계를 수행하는데, 이때는 두 번째 행의 두 번째 열을 고려합니다. 먼저 그 행을 $a_{22}$로 나눈 다음, (두 번째 행을 제외한) 두 번째 열 전체가 $0$이 되도록 다른 행에서 그 행을 뺍니다.
이 과정을 행렬 $A$의 모든 열에 대해 계속합니다. $n = m$이면, $A$는 단위 행렬이 됩니다.
피벗 원소 찾기¶
앞서 설명한 방식에서는 많은 세부 사항이 생략되었습니다. $i$번째 단계에서 $a_{ii}$가 0이라면, 설명한 방법을 바로 적용할 수 없습니다. 대신 먼저 select a pivoting row해야 합니다. 즉, 행렬에서 $i$번째 열이 0이 아닌 행을 하나 찾은 다음, 두 행을 서로 바꿉니다.
여기서는 열이 아니라 행을 교환한다는 점에 유의하십시오. 열을 교환하면 해를 찾은 뒤 원래 위치로 다시 교환하는 것을 잊지 말아야 하기 때문입니다. 따라서 행을 교환하는 편이 훨씬 쉽습니다.
많은 구현에서는 $a_{ii} \neq 0$일 때에도 $a_{ji}$의 절댓값이 가장 큰 행을 피벗 행으로 선택하는 것과 같은 휴리스틱을 사용하여 $i$번째 행을 어떤 피벗 행과 교환하는 것을 볼 수 있습니다. 이 휴리스틱은 이후 단계에서 행렬 값의 범위를 줄이기 위해 사용됩니다. 이 휴리스틱이 없으면 크기가 $20$ 정도인 행렬에서도 오차가 너무 커져 C++의 부동소수점 자료형에서 오버플로가 발생할 수 있습니다.
퇴화된 경우¶
$m = n$이고 시스템이 비퇴화인 경우(즉, 행렬식이 0이 아니고 해가 유일한 경우), 위에서 설명한 알고리즘은 $A$를 단위 행렬로 변환합니다.
이제 $n$과 $m$이 반드시 같지는 않고 시스템이 퇴화될 수도 있는 general case를 살펴봅시다. 이러한 경우에는 $i$번째 단계에서 피벗 원소를 찾지 못할 수 있습니다. 이는 $i$번째 열에서 현재 행부터 아래의 모든 원소가 0이라는 뜻입니다. 이 경우 변수 $x_i$가 가질 수 있는 값이 없거나(즉, SLAE의 해가 없음), $x_i$가 독립 변수여서 임의의 값을 가질 수 있습니다. 가우스-조르당을 구현할 때는 이후의 변수들에 대해 작업을 계속하고 $i$번째 열은 그냥 건너뛰어야 합니다(이는 행렬에서 $i$번째 열을 제거하는 것과 같습니다).
따라서 이 과정에서 일부 변수가 독립 변수로 판명될 수 있습니다. 변수의 개수 $m$이 방정식의 개수 $n$보다 크면 적어도 $m - n$개의 독립 변수가 발견됩니다.
일반적으로 독립 변수를 하나라도 찾으면 그 변수는 임의의 값을 가질 수 있으며, 나머지 (종속) 변수들은 그 변수를 통해 표현됩니다. 이는 실수체에서 작업할 때 시스템이 잠재적으로 무한히 많은 해를 가진다는 뜻입니다. 그러나 독립 변수가 있을 때 SLAE에 해가 전혀 없을 수도 있다는 점을 기억하십시오. 이는 처리되지 않고 남은 방정식 중 적어도 하나의 상수항이 0이 아닐 때 발생합니다. 모든 독립 변수에 0을 대입하고 나머지 변수를 계산한 다음, 원래 SLAE에 대입하여 만족하는지 확인하면 이를 검사할 수 있습니다.
구현¶
다음은 가우스-조르당의 구현입니다. 피벗 행은 현재 열에서 최댓값을 선택하는 휴리스틱으로 고릅니다.
함수 gauss의 입력은 시스템 행렬 $a$입니다. 이 행렬의 마지막 열은 벡터 $b$입니다.
이 함수는 시스템의 해의 개수 $(0, 1,\textrm{or } \infty)$를 반환합니다. 해가 하나 이상 존재하면 그 해는 벡터 $ans$에 담겨 반환됩니다.
const double EPS = 1e-9;
const int INF = 2; // it doesn't actually have to be infinity or a big number
int gauss (vector < vector<double> > a, vector<double> & ans) {
int n = (int) a.size();
int m = (int) a[0].size() - 1;
vector<int> where (m, -1);
for (int col=0, row=0; col<m && row<n; ++col) {
int sel = row;
for (int i=row; i<n; ++i)
if (abs (a[i][col]) > abs (a[sel][col]))
sel = i;
if (abs (a[sel][col]) < EPS)
continue;
for (int i=col; i<=m; ++i)
swap (a[sel][i], a[row][i]);
where[col] = row;
for (int i=0; i<n; ++i)
if (i != row) {
double c = a[i][col] / a[row][col];
for (int j=col; j<=m; ++j)
a[i][j] -= a[row][j] * c;
}
++row;
}
ans.assign (m, 0);
for (int i=0; i<m; ++i)
if (where[i] != -1)
ans[i] = a[where[i]][m] / a[where[i]][i];
for (int i=0; i<n; ++i) {
double sum = 0;
for (int j=0; j<m; ++j)
sum += ans[j] * a[i][j];
if (abs (sum - a[i][m]) > EPS)
return 0;
}
for (int i=0; i<m; ++i)
if (where[i] == -1)
return INF;
return 1;
}
구현 참고 사항:
- 이 함수는 두 개의 포인터, 즉 현재 열 $col$과 현재 행 $row$를 사용합니다.
- 각 변수 $x_i$에 대해 값 $where(i)$는 이 열이 0이 아닌 행입니다. 일부 변수가 독립 변수일 수 있으므로 이 벡터가 필요합니다.
- 이 구현에서는 위에서 설명한 것과 달리 현재 $i$번째 행을 $a_{ii}$로 나누지 않으므로, 마지막에 행렬이 단위 행렬이 되지 않습니다(다만 $i$번째 행을 나누면 오차를 줄이는 데 도움이 되는 것으로 보입니다).
- 해를 찾은 후에는 연립방정식에 해가 적어도 하나 존재하는지 확인하기 위해 그 해를 다시 행렬에 대입합니다. 검증에 성공하면 함수는 독립 변수가 적어도 하나 있는지에 따라 1 또는 $\inf$를 반환합니다.
복잡도¶
이제 이 알고리즘의 복잡도를 추정해 봅시다. 알고리즘은 $m$개의 단계로 이루어지며, 각 단계에서는 다음을 수행합니다.
- 피벗 행을 탐색하고 행의 순서를 바꿉니다. 앞서 언급한 휴리스틱을 사용하면 이 작업에 $O(n + m)$이 걸립니다.
- 현재 열에서 피벗 원소를 찾았다면 이 방정식을 다른 모든 방정식에 더해야 하며, 여기에는 $O(nm)$의 시간이 걸립니다.
따라서 알고리즘의 최종 복잡도는 $O(\min (n, m) . nm)$입니다. $n = m$인 경우 복잡도는 간단히 $O(n^3)$입니다.
SLAE가 실수 위에서가 아니라 모듈로 2에서 주어진 경우에는 연립방정식을 훨씬 더 빠르게 풀 수 있으며, 이는 아래에서 설명합니다.
알고리즘의 가속¶
알고리즘을 전진 단계와 후진 단계의 두 단계로 나누면 앞의 구현을 두 배 빠르게 만들 수 있습니다.
- 전진 단계: 앞의 구현과 비슷하지만, 현재 행을 그 뒤에 있는 행들에만 더합니다. 그 결과 대각 행렬 대신 삼각 행렬을 얻습니다.
- 후진 단계: 행렬이 삼각 행렬이면 먼저 마지막 변수의 값을 계산합니다. 그런 다음 이 값을 대입하여 그다음 변수의 값을 구합니다. 그런 다음 이 두 값을 대입하여 그다음 변수들을 구합니다...
후진 단계는 $O(nm)$만 걸리므로 전진 단계보다 훨씬 빠릅니다. 전진 단계에서는 연산 횟수를 절반으로 줄이므로 구현의 실행 시간이 줄어듭니다.
모듈러 SLAE 풀기¶
어떤 모듈러스에서 SLAE를 풀 때에도 앞에서 설명한 알고리즘을 그대로 사용할 수 있습니다. 그러나 모듈러스가 2인 경우에는 비트 연산과 C++ bitset 자료형을 사용하여 가우스-조르당 소거법을 훨씬 더 효율적으로 수행할 수 있습니다.
int gauss (vector < bitset<N> > a, int n, int m, bitset<N> & ans) {
vector<int> where (m, -1);
for (int col=0, row=0; col<m && row<n; ++col) {
for (int i=row; i<n; ++i)
if (a[i][col]) {
swap (a[i], a[row]);
break;
}
if (! a[row][col])
continue;
where[col] = row;
for (int i=0; i<n; ++i)
if (i != row && a[i][col])
a[i] ^= a[row];
++row;
}
// The rest of implementation is the same as above
}
비트 압축을 사용하므로 구현이 더 짧을 뿐만 아니라 32배 더 빠릅니다.
피벗 행을 고르는 여러 휴리스틱에 대한 간단한 참고 사항¶
어떤 휴리스틱을 사용해야 하는지에 대한 일반적인 규칙은 없습니다.
앞의 구현에서 사용한 휴리스틱은 실제로 꽤 잘 동작합니다. 또한 이 휴리스틱은 (현재 행과 현재 열부터 시작하는) 부분 행렬 전체의 모든 원소 중에서 피벗 행을 찾는 "완전 피벗팅"과 거의 같은 답을 내는 것으로 밝혀졌습니다.
다만 두 휴리스틱 모두 원래 방정식들이 얼마나 스케일링되었는지에 의존한다는 점에 유의하십시오. 예를 들어 방정식 중 하나에 $10^6$을 곱했다면, 이 방정식은 첫 번째 단계에서 거의 확실하게 피벗으로 선택됩니다. 이는 다소 이상해 보이므로, implicit pivoting이라고 부르는 더 복잡한 휴리스틱으로 바꾸는 것이 타당해 보입니다.
암묵적 피벗팅은 두 행이 모두 정규화되어 최대 원소가 1이 된 것처럼 원소를 비교합니다. 이 기법을 구현하려면 각 행의 최댓값을 유지해야 합니다(또는 최댓값이 1이 되도록 각 행을 유지할 수도 있지만, 이렇게 하면 누적 오차가 커질 수 있습니다).
해 개선하기¶
여러 가지 휴리스틱을 사용하더라도, 가우스-조르당 알고리즘은 크기가 $50 - 100$에 불과한 특수한 행렬에서도 여전히 큰 오차를 낼 수 있습니다.
따라서 가우스-조르당으로 얻은 해는 때때로 간단한 수치적 방법, 예를 들어 단순 반복법을 적용하여 개선해야 합니다.
따라서 풀이는 두 단계로 이루어집니다. 먼저 가우스-조르당 알고리즘을 적용하고, 그다음 첫 번째 단계에서 얻은 해를 초기 해로 삼아 수치적 방법을 적용합니다.