The Story of Joon
Linear Algebra in Problem Solving (1) 본문
Linear Algebra in Problem Solving (1) (현 포스트)
Linear Algebra in Problem Solving (2)
Linear Algebra in Problem Solving (3)
선형대수학(linear algebra)은 벡터 공간(vector space)과 선형 변환(linear transformation)에 대해서 다루는 학문으로, 여러 분야에서 많은 응용이 되고 있고 수학에서의 중요성 또한 매우 높다. 컴퓨터 과학도 예외는 아닌데, 이런 이유로 많은 학교에서 컴퓨터 과학/공학과에서 기초 선형대수학을 필수 과목으로 지정하기도 한다.
이 포스트 시리즈에서는 PS에서 등장하는 선형대수학 관련 문제들을 해결하기 위한 여러 가지 선형대수학 알고리즘을 다룰 예정이다. 일반적인 PS에서 선형대수학이 필요한 경우는 드물지만, BOJ에 관련 문제들을 종종 찾아볼 수 있고 난이도가 높은 문제 셋에서는 상당한 수준의 선형대수학 배경지식을 요구하는 문제를 어렵지 않게 만날 수 있다. 그러한 문제에 접근할 수 있는 알고리즘적 토대를 제공하는 것이 이 시리즈의 목적이다.
아쉽게도 이 시리즈에서는 기초적인 선형대수학 개념을 다루지는 않는다. 선형대수학을 토대부터 쌓는 것은 굉장히 시간이 많이 걸리기 때문이다. 대신 필요할 경우 위키피디아 등 잘 정리된 레퍼런스로 링크할 예정이다. 이 글을 무리없이 읽기 위해서는 대학 1~2학년 수준의 선형대수학 지식을 갖추는 것이 권장된다. 추후에 이 블로그에서 기초 선형대수학을 다룰 생각이 없지는 않지만, 당분간은 계획이 없다.
본 시리즈의 첫 번째인 이 포스트에서는 행렬 곱셈, 가우스 소거법, 역행렬, 행렬식 등 기초적인 값을 구하는 방법에 대해 다룬다.
언어는 C++로 진행하며, 편의상 아래와 같은 클래스를 바탕으로 설명한다. 설명의 편의를 위해 간략하게 만든 코드이며, 필자도 문제를 풀 때 이 클래스를 쓰지는 않는다. 실제 구현할 때는 본인이 원하는 형태로 구현하면 된다.
template<class T>
class matrix {
public:
int n, m;
vector<vector<T>> entry;
matrix(int n_, int m_) : n(n_), m(m_), entry(n_, vector<T>(m_)) {}
matrix(int n_) : matrix(n_, n_) {}
T& operator[](const pair<int, int>& p) {
return entry[p.first][p.second];
}
const T& operator[](const pair<int, int>& p) const {
return entry[p.first][p.second];
}
};
행렬의 원소
기초 선형대수학을 공부해본 독자라면 행렬의 원소가 실수(\(\mathbb{R}\))고 벡터공간 또한 \(\mathbb{R}^n\)인 경우가 가장 익숙할 것이다. 하지만 대부분의 선형대수학 이론은 실수보다 더욱 일반적인 체(field)라는 대수적 구조에서도 전개가 가능하다. 심지어 체에서 나눗셈이 없어진 환(ring)이라는 더욱 일반화된 대수적 구조에서도 전개가 가능한 이론도 많지만, 이쪽은 선형대수학이 아니라 모듈(module) 이론이라고 부른다.
체는 간단하게 말해서 덧셈, 곱셈 연산이 있고, 각각에 대한 항등원과 역원(단, 덧셈에 대한 항등원은 곱셈에 대한 역원의 요구사항에서 배제된다)이 있으며 교환법칙, 결합법칙, 그리고 분배법칙이 성립하는 대수적 구조이다. 대표적으로 우리에게 가장 익숙한 실수 집합 \(\mathbb{R}\)과 복소수 집합 \(\mathbb{C}\)가 있으며, 정수를 소수 \(p\)로 나눈 나머지들의 집합인 \(\mathbb{F}_p\)도 있다. 특히 \(p=2\)인 경우는 우리에게 익숙한 비트의 집합이다.
실수 집합과 복소수 집합의 경우 실생활에서 보기 쉽고 이론을 전개하기에 편한 점이 있지만, 컴퓨터에서 구현할 경우 사칙연산이 상당히 많은 선형대수 알고리즘 특성상 실수 오차로 인해 수치적으로 불안정해지고 정확한 결과를 기대하기 어려운 경우가 많다. 따라서 PS에서 보기 가장 흔한 체는 \(\mathbb{F}_p\)이다. 예를 들어 선형대수학으로 풀 수 있는 문제가 결과를 998244353으로 나눈 나머지를 출력하라고 한다면, \(p=998244353\)인 체에서 문제를 풀면 된다.
요약하자면, 행렬의 원소가 반드시 실수일 필요가 없음을 기억하고 있어야 한다.
행렬 곱셈
행렬 곱셈은 가장 기본적인 행렬 연산으로 아래와 같은 코드로 쉽게 해결할 수 있다.
template<class T>
matrix<T> matmul(const matrix<T>& A, const matrix<T>& B) {
assert(A.m == B.n);
matrix<T> C(A.n, B.m);
for (int i = 0; i < C.n; i++) {
for (int k = 0; k < A.m; k++) {
for (int j = 0; j < C.m; j++) {
C[{i, j}] += A[{i, k}] * B[{k, j}];
}
}
}
return C;
}여기서 다중 반복문을 돌리는 순서에 따라 성능이 달라질 수 있는데, 이 코드처럼 i-k-j 순서로 돌리는 것이 i-j-k 순서로 돌리는 것보다 더 캐시 친화적이다.
이렇게 구현하면 \(n\times m\) 행렬과 \(m\times k\) 행렬을 곱할 때 \(O(nmk)\) 시간이 걸린다. 이것보다 더 빠른 행렬 곱셈 중 그나마 실제로 구현해 볼만한 것은 Strassen algorithm인데, PS에서 활용할 범위 내에서는 성능 향상이 그리 크지 않을 뿐더러 작은 행렬에 대해서는 오히려 나이브한 곱셈이 더 빠르기 때문에 그다지 권장하지 않는다. 그래도 관심이 있는 독자는 이 포스트를 추천한다.
만약 \(\mathbb{F}_2\) 위의 행렬, 즉 행렬의 원소가 비트라면 bitset을 활용해 좀더 효율적으로 구현할 수 있다. 단 POPCNT 연산이 하드웨어의 지원을 받는 경우에 한한다.
template<size_t M, size_t K>
vector<bitset<M>> matmul(const vector<bitset<M>>& A, const vector<bitset<K>>& B) {
int n = (int)A.size(), m = (int)B.size();
assert(m == M);
vector<bitset<K>> C(n);
vector<bitset<M>> Bt(K);
for (int i = 0; i < K; i++) {
for (int j = 0; j < m; j++) {
Bt[i][j] = B[j][i];
}
}
for (int i = 0; i < n; i++) {
for (int j = 0; j < K; j++) {
C[i][j] = (A[i] & Bt[j]).count() & 1;
}
}
return C;
}이때 시간복잡도는 bitset의 word size를 \(w\)라고 했을 때 \(O(nmk/w)\)로, 일반적으로 \(w=64\)임을 생각해 본다면 꽤 큰 성능 향상을 기대해볼 수 있다.
가우스 소거법
가우스 소거법(Gaussian elimination)은 앞으로 나올 대부분의 선형대수 알고리즘의 기본 토대이므로 반드시 능숙하게 다룰 수 있어야 한다.
구현은 우리에게 익숙한 가우스 소거법 과정을 그대로 시뮬레이션하면 된다. 가우스 소거법 과정을 복습해보자.
- 1열(column)에서 시작해서 leading one의 위치를 찾기 시작한다.
- 이전 leading one이 위치한 다음 행(row)부터 차례대로 보기 시작해서 leading one으로 선정할 행을 찾는다. 즉 해당 위치에 0이 아닌 원소가 있어야 한다.
- 그러한 행이 없다면 그 열은 leading one이 없는 열이다. 다음 열로 넘어가서 2번 과정으로 돌아간다.
- 그러한 행을 찾았으면 이전 leading one 바로 다음 행에 위치하도록 위치를 교환한다. (Row switching)
- 적절한 상수를 곱해 leading one 위치에 1이 오도록 한다. (Row multiplication)
- 해당 행에 적절한 상수를 곱하고 그 다음 행들에 더해 leading one의 아래쪽으로는 모두 0이 되도록 한다. (Row addition)
참고로 row switching, row multiplication, row addition을 한 데 묶어 elementary row operation이라고 하고, 꽤 중요한 개념이니 알아두면 좋다.
이제 이를 직접 구현하면 된다. 우선 편의를 위해 행렬 클래스에 elementary row operation을 수행하는 함수를 만든다. 여기서 row switching 연산은 상수 시간에 효율적으로 구현할 수 있는데, 직접 원소들을 바꿔줄 필요 없이 레퍼런스를 바꿔주면 된다.
template<class T>
class matrix {
public:
int n, m;
vector<vector<T>> entry;
matrix(int n_, int m_) : n(n_), m(m_), entry(n_, vector<T>(m_)) {}
matrix(int n_) : matrix(n_, n_) {}
T& operator[](const pair<int, int>& p) {
return entry[p.first][p.second];
}
const T& operator[](const pair<int, int>& p) const {
return entry[p.first][p.second];
}
void rswap(int i, int j) {
swap(entry[i], entry[j]);
}
void rmul(int i, T x) {
for (int j = 0; j < m; j++) {
entry[i][j] *= x;
}
}
void radd(int i1, int i2, T x) {
for (int j = 0; j < m; j++) {
entry[i1][j] += entry[i2][j] * x;
}
}
};이를 활용해서 가우스 소거법을 구현하면 된다.
template<class T>
void gaussian_elimination(matrix<T>& A) {
for (int j = 0, r = 0; j < A.m && r < A.n; j++) {
for (int i = r; i < A.n; i++) {
if (A[{i, j}] != 0) {
A.rswap(r, i);
break;
}
}
if (A[{r, j}] != 0) {
A.rmul(r, 1 / A[{r, j}]);
for (int i = r + 1; i < A.n; i++) {
A.radd(i, r, -A[{i, j}]);
}
r++;
}
}
}\(n\times m\) 행렬의 가우스 소거법을 계산하는 시간복잡도는 \(O(nm\cdot\min(n, m))\). 행렬 곱셈에서와 마찬가지로, \(\mathbb{F}_2\) 위라면 bitset을 활용할 수 있다.
template<size_t M>
void gaussian_elimination(vector<bitset<M>>& A) {
int n = (int)A.size();
for (int j = 0, r = 0; j < M && r < n; j++) {
for (int i = r; i < n; i++) {
if (A[i][j]) {
swap(A[r], A[i]);
break;
}
}
if (A[r][j]) {
for (int i = r + 1; i < n; i++) {
if (A[i][j]) {
A[i] ^= A[r];
}
}
r++;
}
}
}이때 시간복잡도는 \(O(nm\cdot\min(n, m)/w)\).
행렬의 rank 구하기
가우스 소거법의 부산물로 행렬의 rank도 구할 수 있다. 가우스 소거법이 완료된 후 leading one의 개수가 곧 행렬의 rank이므로, 위 코드에서 변수 r에 행렬의 rank가 저장된다.
template<class T>
int rank(matrix<T>& A) {
int r = 0;
for (int j = 0; j < A.m && r < A.n; j++) {
for (int i = r; i < A.n; i++) {
if (A[{i, j}] != 0) {
A.rswap(r, i);
break;
}
}
if (A[{r, j}] != 0) {
A.rmul(r, 1 / A[{r, j}]);
for (int i = r + 1; i < A.n; i++) {
A.radd(i, r, -A[{i, j}]);
}
r++;
}
}
return r;
}행렬식
\(n\times n\) 정사각 행렬 \(A\)의 행렬식은 아래와 같은 식의 값으로 알려져 있다.
\[\det(A)=\sum_{\sigma\in S_n}\left(\operatorname{sgn}(\sigma)\prod_{i=1}^n a_{i,\sigma_i}\right)\]
여기서 \(S_n\)은 \(\{1,\cdots,n\}\)의 모든 순열의 집합이고, \(\operatorname{sgn}\)은 순열의 홀짝성을 나타내는 함수이다. 순열의 홀짝성이란 간단히 말해 두 원소의 swap만으로 순열을 정렬할 때 필요한 swap 횟수의 홀짝성이다. (어떤 과정을 거치든 홀짝성은 같다.) 이 횟수가 짝수면 \(\operatorname{sgn}\)의 값이 1이고, 그렇지 않을시 -1이다.
행렬식을 계산할 때 나이브하게 구하는 것 외에 cofactor expansion이라는 방식도 있지만, 모두 지수 시간이 걸리는 비효율적인 알고리즘이다. 하지만 위에서 구현한 가우스 소거법을 이용하면 \(O(n^3)\)에 행렬식을 구할 수 있다. 가우스 소거법이 완료되면 주대각선 밑으로는 모두 0인 상삼각행렬(upper triangular matrix)이 되는데, 이런 행렬의 행렬식은 주대각선의 원소를 모두 곱하면 되므로 쉽게 구할 수 있기 때문이다.
따라서 각각의 elementary row operation이 행렬식을 어떻게 변화시키는지 보면 되는데, 매우 간단하다.
- Row switching: 행렬식의 부호만 뒤집힌다.
- Row multiplication: 행렬식에도 곱한 상수만큼 곱해진다.
- Row addition: 변화가 없다.
따라서 아래와 같이 구현하면 된다.
template<class T>
T det(matrix<T>& A) {
assert(A.n == A.m);
T ret = 1;
for (int j = 0; j < A.n; j++) {
for (int i = j; i < A.n; i++) {
if (A[{i, j}] != 0) {
if (i > j) {
A.rswap(j, i);
ret = -ret;
}
break;
}
}
if (A[{j, j}] == 0) {
return 0;
}
ret *= A[{j, j}];
A.rmul(j, 1 / A[{j, j}]);
for (int i = j + 1; i < A.n; i++) {
A.radd(i, j, -A[{i, j}]);
}
}
return ret;
}Reduced row echelon form (RREF)
가우스 소거법의 결과로 만들어지는 행렬 꼴을 row echelon form이라고 하는데, reduced row echelon form은 그것보다 좀더 처리가 된 행렬이다.
단순 row echelon form이 각각의 leading one 아래로는 모두 0인 행렬이었다면, RREF는 단순히 그 위까지도 0인 행렬을 의미한다. 따라서 구현도 반복문 하나만 수정해주면 된다.
template<class T>
void rref(matrix<T>& A) {
for (int j = 0, r = 0; j < A.m && r < A.n; j++) {
for (int i = r; i < A.n; i++) {
if (A[{i, j}] != 0) {
A.rswap(r, i);
break;
}
}
if (A[{r, j}] != 0) {
A.rmul(r, 1 / A[{r, j}]);
for (int i = 0; i < A.n; i++) {
if (i != r) {
A.radd(i, r, -A[{i, j}]);
}
}
r++;
}
}
}선형 연립 방정식 풀기
선형 연립 방정식은 \(n\times m\) 행렬 \(A\)와 크기 \(n\)인 벡터 \(\mathbf{b}\)가 있을 때, 아래 식을 만족하는 크기 \(m\)인 벡터 \(\mathbf{x}\)를 구하는 문제이다.
\[A\mathbf{x}=\mathbf{b}\]
가장 간단한 형태는 \(n=m\)이고 \(A\)가 역행렬이 존재하는 경우이다. 그래야 아래와 같이 해가 유일하기 때문이다.
\[\mathbf{x}=A^{-1}\mathbf{b}\]
개념상으로는 이렇지만 실제로 역행렬을 구해서 이 문제를 해결하는 것은 권장되지 않는데, 이유는 여러 가지가 있다. 첫째로 수치해석학적인 관점으로 볼 때 역행렬을 거쳐서 해를 구하게 되면 수치적으로 불안정하고, 두 번째로는 역행렬을 구하는게 이 문제를 푸는 것보다 어렵기 때문이다.
권장되는 방법은 \(A\)에 \(\mathbf{b}\)를 붙여서 augmented matrix를 만든 후 RREF 형태를 만드는 것이다.
\[\begin{bmatrix} A & |\mathbf{b}\end{bmatrix}\]
\(A\)의 역행렬이 존재한다는 조건으로 인해 RREF 형태가 되면 원래 \(A\)였던 부분은 반드시 단위행렬이 되고, \(\mathbf{b}\) 부분에는 우리가 구하고 싶었던 \(\mathbf{x}\)가 남게 된다. 이는 우리가 선형 연립방정식을 손으로 푸는 과정과 일치한다.
template<class T>
vector<T> solve_linear_equation(const matrix<T>& A, const vector<T>& b) {
int n = A.n;
assert(n == A.m);
assert(n == (int)b.size());
matrix<T> B(n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
B[{i, j}] = A[{i, j}];
}
B[{i, n}] = b[i];
}
rref(B);
vector<T> x(n);
for (int i = 0; i < n; i++) {
assert((B[{i, i}] == 1));
x[i] = B[{i, n}];
}
return x;
}\(A\)의 역행렬이 존재하지 않거나 \(n\neq m\)인 경우에도 마찬가지 방법을 사용할 수 있다. 이 경우 RREF 변환이 끝난 후 어떤 결과가 나오느냐에 따라 해가 존재하는지, 존재하다면 유일한지 혹은 무한한지 등을 판별할 수 있다. 독자들은 선형대수학을 이미 어느 정도 알고 있으므로 이는 쉬운 연습문제로 남긴다.
응용: XOR이 특정 수가 되는 부분집합 구하기
\(2^{w}\)보다 작은 음이 아닌 정수 \(n\)개와 \(b\)가 주어질 때, \(n\)개의 수 중 XOR을 했을 때 \(b\)가 되는 수의 부분집합을 찾는 문제를 생각해보자. 만약 XOR이 아니라 합이었다면 그 유명한 subset sum 계열의 문제가 되고, NP-complete이다. 하지만 XOR의 경우 선형대수학과 bitset을 활용해 \(O(wn\cdot\min(w, n)/w)\cong O(wn)\)에 해결할 수 있다.
핵심은 각 수를 \(\mathbb{F}_2\)에서의 \(w\)차원 벡터로 취급하는 것이다. 그러면 XOR은 이 벡터공간에서는 합이 되기 때문에, 결국 합이 특정 벡터가 되는 벡터 집합을 찾는 문제가 된다. \(\mathbb{F}_2\)에서는 선형 결합이 곧 어떤 부분집합의 원소들의 합이기 때문에, 주어진 \(n\)개의 벡터의 선형결합이 \(b\)가 될 수 있는지, 될 수 있다면 그 선형결합을 어떻게 구성해야 하는지 찾는 문제이다. 즉 선형 연립방정식을 푸는 것과 같다.
이렇듯 음이 아닌 정수를 \(\mathbb{F}_2\) 위에서의 벡터, XOR을 벡터의 합으로 해석해서 선형대수학적으로 접근하는 문제들을 흔히 만날 수 있기 때문에 알아두면 좋다. 선형대수 알고리즘을 직접 사용하지 않더라도 문제에 접근하는데 많은 도움을 준다.
역행렬
역행렬은 문제에서 대놓고 요구하지 않는 이상 직접 구할 일이 흔하지는 않지만, 선형 연립 방정식의 연장선으로서 알아두어서 나쁠 것 없을 것 같아서 소개한다.
\(n\times n\) 행렬 \(A\)의 역행렬을 구하기 위해서는 아래 식을 만족하는 행렬 \(X\)를 구하면 된다.
\[AX=I\]
\(I\)는 \(n\times n\) 단위행렬로, 주대각선은 모두 1이고 나머지는 0인 행렬이다.
역행렬을 구하는 방법으로 Cramer's rule 같은 것을 떠올릴 수 있지만, 행렬식을 구할 때와 마찬가지로 지수 시간이 걸리므로 실제로 사용할 수는 없다.
이 문제를 열 단위로 쪼개서 생각해보자. \(X\)의 \(i\)번째 열을 구하기 위해서는 아래 식을 만족하는 벡터 \(\mathbf{x}_i\)를 구하면 된다.
\[A\mathbf{x}_i=\mathbf{\delta}_i\]
여기서 \(\mathbf{\delta}_i\)는 \(i\)번째 원소만 1이고 나머지는 0인 벡터이다. 이 문제를 선형 연립 방정식을 풀 때와 동일한 환경이다. 따라서 각 \(\mathbf{\delta}_i\)에 대해 augmented matrix를 만든 후 가우스 소거법을 \(n\)번 실행하면 \(X\)의 각 열을 구할 수 있지만, \(A\) 부분이 소거되는 과정은 완전히 동일하므로 굳이 이렇게 할 필요 없다. 아예 \(I\)를 통으로 augment하면 된다.
\[\begin{bmatrix}A&|&I\end{bmatrix}\]
이후 RREF 형태로 만들면 \(A\) 부분은 단위행렬이 되고 \(I\) 부분은 \(A\)의 역행렬이 된다.
template<class T>
matrix<T> matinv(const matrix<T>& A) {
int n = A.n;
assert(n == A.m);
matrix<T> B(n, n + n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
B[{i, j}] = A[{i, j}];
B[{i, n + i}] = 1;
}
}
rref(B);
matrix<T> X(n);
for (int i = 0; i < n; i++) {
assert((B[{i, i}] == 1));
for (int j = 0; j < n; j++) {
X[{i, j}] = B[{i, n + j}];
}
}
return X;
}마치며
이 포스트에서는 선형대수학의 가장 기본적인 알고리즘을 알아보았다. 다음 포스트에서는 좀더 심화된 선형대수 알고리즘(특성다항식, 일차 다항식 행렬, 종결식 등)과 여러 응용(그래프 이론, 매트로이드, 조합론 등)에 대해서 다룰 예정이다.
연습문제
난이도 순이 아니다.
BOJ 13712: 키르히호프의 법칙 (회로 이론에서의 키르히호프 법칙에 대한 사전지식이 필요하다.)
BOJ 15757: 일반 그래프 매칭 (Tutte matrix에 대한 사전지식이 필요하다.)
'Computer Science > 알고리즘' 카테고리의 다른 글
| Linear-time construction of the sieve of Eratosthenes (0) | 2022.09.11 |
|---|---|
| Linear Algebra in Problem Solving (2) (2) | 2022.09.02 |
| 2015 ACM-ICPC 한국 예선 F - 파일 합치기 (2) | 2020.05.02 |
| Link/Cut Tree (2) (7) | 2017.08.21 |
| Fast Fourier Transform (8) | 2017.08.15 |
