이상엽 선형대수학 정주행 2회차 - 행렬식과 역행렬

1. 소행렬과 행렬식

1) 소행렬(Submatrix): 특정 행과 열을 제거했을 때, 만들어지는 행렬

행렬식을 알기 위해서는 소행렬이 뭔지 알아야 한다. 소행렬은 어떤 행렬에서 선택된 행과 열을 제외하고 만들어지는 행렬을 말한다. 예를 들어 행렬 A의 소행렬 M(i, j)은 아래 이미지와 같이 나오는 것이다. 참고로 여기서 얻는 행렬식을 소행렬식이라고 부른다.

소행렬 M(1, 1)이라고 할 수 있다.

 

2) 행렬식(Determinant): 행렬을 하나의 값에 대응시키는 함수

행렬식은 어떤 정사각행렬을 대표하는 값을 구하기 위한 하나의 함수라고 생각하면 된다. 표기는 det(A) = |A| 와 같이 표기하며 2x2까지는 쉽게 |A|의 값을 구할 수 있지만 3x3부터 번거로워진다.

 

0x0: 행렬의 길이가 0인 경우 즉, 행렬의 성분이 없는 이런 경우는 0으로 정의되어 있다.

1x1: 행렬의 원소가 1개이면 이 행렬을 대표할 수 있는 값 역시 하나 밖에 없다. 따라서 해당 성분의 값 자체가 |A| 이다.

 

2x2부터는 다음과 같은 과정으로 행렬식이 이루어진다. 읽어도 아마 이해안될거다. 한 번 읽어나보고 이미지와 함께 아래의 설명과 같이 읽으면 그나마 이해할 수 있다. 어쨋든 행렬 A가 있다고 가정하면

1. 1번째 행의 1번째 열부터 시작하며 3~4의 과정을 재귀적으로 반복하면 된다. (사실 굳이 1번째 행이 아니어도 됨)
2. 앞으로 현재 행을 i, 현재 열을 j라고 가정하고 진행할 것이다.
3. A(i, j)인 성분과 해당 성분이 속해있는 행과 열을 제외한 소행렬식(소행렬의 행렬식)을 곱한다. (M(i, j)로 표기)
4. 이제 (-1)^(i + j) x M(i, j)의 값을 더한다. 이걸 여인수라고 하며 나중에 이 값의 총합을 구할 것이다.
5. 현재 열이 다음 열로 넘어가고 3번으로 돌아간다. 다음 열로 넘어갈 것이 없으면 연산을 끝낸다.
6. 최종적으로 4번 과정을 통해 전부 더해진 값이 행렬식의 값이다.

꽤 복잡하게 느껴질 것 같기도 한데 그렇다면 아래의 그림을 보자.

3x3 행렬식의 과정

지금 보고있는 이미지는 3x3 행렬식의 과정이다. 그런데 뜬금없이 2x2의 행렬들로 쪼개지는 것을 볼 수 있을 것이다. 사실 저게 뭐냐면 2x2 행렬식을 표현한 것이다. 이게 위에서 서술한 2의 과정이라고 보면 된다. 3x3 행렬식을 위해 2x2 행렬식을 진행시키는 것이다. 그럼 2x2 행렬식의 과정은 어떻게 될까? 위에서 서술한 것을 똑같이 하면 된다. 그럼 1x1 행렬식 즉, 하나의 값이 나올텐데 그게 곧 2x2 행렬식의 과정에 사용되는 수인 것이다.

 

그럼 4x4 행렬식의 과정은 어떻게 될까? 사람이 직접 계산하기에는 계산 과정이 워낙 실수가 나올 여지가 많고 번거로워서 그렇지 연산을 하는 순서 즉, 개념상으로는 그렇게까지 어렵지 않을 것이다. a11~a44까지의 성분이 4x4의 형태로 배열되어 있다면 위의 그림과 똑같은 원리로 a11~a14와 a11~a14에 해당하는 소행렬식을 이용하면 된다. 둘을 곱해서 3번에서 서술한 것처럼 더하고 빼면 되는 것이다.

 

사실 이 알고리즘의 시간 복잡도는 무려 n! 이나 된다....ㄷㄷ 아무리 생각해도 코드로 구현하기에는 시간 복잡도가 너무 치명적인 것 같아서 다른 알고리즘이 없는지 찾아보았다. 그런데 막상 알고리즘을 들어보니까 뜬금없이 기하학적인 내용들이 우후죽순 튀어나왔다. 갑자기 직교가 어쩌고 투영이 어쩌고... 뭔가 이상해서 이상엽Math 강의 목록을 다시 살펴봤는데 행렬 바로 다음에 벡터가 나오는 것을 알 수 있었다. 아마 더 효율적인 알고리즘으로 행렬식을 구하는 방법은 그 부분을 배우고 나서 배울 수 있는 듯 하다. 일단 이번 포스팅에서는 지금 배운 알고리즘으로 구현을 해보도록 하겠다.

 

 

2. 수반행렬과 역행렬

1) 수반행렬(Adjoint of Matrix): 여인수로 이루어진 행렬의 전치행렬

역행렬을 알기 위해서 수반행렬이 뭔지 알아야 하고 수반행렬을 구하려면 여인수가 뭔지 알아야 한다. 그런데 여인수는 이미 행렬식을 구할 때, 언급한 적이 있다. 행렬식을 구하는 과정 중 4번에서 구한 값이 바로 여인수다. 이제 이 여인수를 이용해서 수반행렬을 구할 것이다. 수반행렬의 정의는 아래의 이미지와 같다.

수반행렬의 정의

여기서 주목해야 할 점은 행렬 A와 수반행렬 adj(A)를 곱하면 모든 대각성분이 행렬식 det(A)인 대각행렬이 나온다는 것이다. 행렬식을 구하는 과정을 이해했다면 왜 이렇게 나오는지 이해하게 될 것이다. 직접 행렬식을 구해보는 것도 지금 이 말을 이해하는 것에 큰 도움이 될 것이다. 어쨋든 지금 아래와 같은 상황이라는 것이다.

 

D = 모든 대각성분이 det(A)인 대각행렬

 

A * adj(A) = D

 

 

2) 역행렬(Inverse Matrix): 어떤 정사각행렬의 곱셈의 역원이 되는 행렬

여기까지 했으면 다왔다. 이제 역행렬을 구할 수 있다. 위에서 D의 모든 대각성분이 det(A)인 대각행렬이라고 했다. 그럼 나머지 성분은 전부 0이니 여기에 det(A)를 상수배로 나누면 어떻게 될까? 단위행렬이 된다. 여기서 역행렬이 어떤 행렬의 곱셈의 역원이 되는 행렬이라고 했다. 이게 무슨 뜻이냐? 어떤 행렬의 곱했는데 그게 단위행렬이 나왔다? 그럼 그건 어떤 행렬의 역행렬이라는 뜻이다. 이제 아래를 보자.

 

A * adj(A) = D

 

D / det(A) = I (단위행렬)

 

A * adj(A) / det(A) = D / det(A)

 

A * adj(A) / det(A) = I (단위행렬)

 

이제 역행렬의 정의에 따라서 adj(A) / det(A)가 역행렬이라는 것을 알아냈다. 즉, 역행렬은 (수반행렬 / 행렬식)이었던 것이다! 참고로 행렬식이 0이 나오는 경우가 있는데 해당 경우는 역행렬을 구하는 것이 불가능하다. 역행렬이 존재하지 않기 때문이다. 그럼 이제 구현하는 코드를 짜볼텐데 기존에 구현했던 것과 같이 함께 올리도록 하겠다.

 

 

3. 번외

강의에서 문제로 몇 개의 공식들을 소개했다. 증명하는 과정까지 있었는데 이 글에서 증명은 생략하겠다.

(A^-1은 A의 역행렬, Ac는 A의 수반행렬)

1. det(A) = 1 / det(A^-1)
2. (A^-1)c = Ac^-1
3. (AB)^-1 = (B^-1)(A^-1)

 

 

4. 구현

Matrix.h

#ifndef __MATRIX_H__
#define __MATRIX_H__

#include <initializer_list>

#define DEFAULT_EPSILON 0.000001f

class Matrix
{
public:

	Matrix(const Matrix& matrix);

	Matrix(const std::initializer_list<std::initializer_list<float>>& matrix);

	~Matrix();

public:

	static Matrix CreateZeroMatrix(size_t rowNum, size_t columnNum);

	static Matrix CreateIdentityMatrix(size_t length);

	size_t GetRowNum() const;
	size_t GetColumnNum() const;

	void SetEpsilon(float epsilon);
	float GetEpsilon() const;

	void GetDiagonalElements(float diagonalElements[]) const;

	Matrix GetTransposeMatrix() const;

	Matrix GetSubmatrix(size_t row, size_t column) const;

	float GetDeterminant() const;

	Matrix GetAdjointMatrix() const;

	Matrix GetInverseMatrix() const;

	bool IsZeroMatrix() const;

	bool IsSymmetricMatrix() const;

	bool IsSkewSymmetricMatrix() const;

	bool IsSquareMatrix() const;

	bool IsDiagonalMatrix() const;

	bool IsIdentityMatrix() const;

public:

	Matrix& operator=(const Matrix& matrix);

	Matrix operator+(const Matrix& matrix) const;
	Matrix operator-(const Matrix& matrix) const;
	Matrix operator*(const Matrix& matrix) const;
	Matrix operator/(const Matrix& matrix) const;

	Matrix operator*(float value) const;
	Matrix operator/(float value) const;
	friend Matrix operator*(float value, const Matrix& matrix);
	friend Matrix operator/(float value, const Matrix& matrix);

	Matrix& operator+=(const Matrix& matrix);
	Matrix& operator-=(const Matrix& matrix);
	Matrix& operator*=(const Matrix& matrix);
	Matrix& operator/=(const Matrix& matrix);

	Matrix& operator*=(float value);
	Matrix& operator/=(float value);

	float* operator[](int i) const;

private:

	size_t rowNum;
	size_t columnNum;

	float** elements;

	float epsilon;
};

#endif

 

Matrix.cpp

#include "Matrix.h"

#include <cmath>
#include <cassert>

Matrix::Matrix(const Matrix& matrix)
	: epsilon(DEFAULT_EPSILON)
	, rowNum(matrix.rowNum)
	, columnNum(matrix.columnNum)
{
	int i;
	int j;

	elements = new float*[rowNum];

	for (i = 0; i < rowNum; ++i)
	{
		elements[i] = new float[columnNum];

		for (j = 0; j < columnNum; ++j)
		{
			elements[i][j] = matrix[i][j];
		}
	}
}

Matrix::Matrix(const std::initializer_list<std::initializer_list<float>>& matrix)
	: epsilon(DEFAULT_EPSILON)
{
	int i = 0;
	int j = 0;

	rowNum = matrix.size();
	columnNum = rowNum <= 0 ? 0 : matrix.begin()->size();
	elements = new float*[rowNum];

	for (std::initializer_list<float> row : matrix)
	{
		elements[i] = new float[columnNum];

		j = 0;

		for (float element : row)
		{
			elements[i][j] = element;
			++j;
		}

		++i;
	}
}

Matrix::~Matrix()
{
	for (int i = 0; i < rowNum; ++i)
	{
		delete[] elements[i];
	}

	delete[] elements;
}

Matrix Matrix::CreateZeroMatrix(size_t rowNum, size_t columnNum)
{
	Matrix matrix = { };

	matrix.rowNum = rowNum;
	matrix.columnNum = columnNum;
	matrix.elements = new float*[rowNum];

	for (int i = 0; i < rowNum; ++i)
	{
		matrix.elements[i] = new float[columnNum] { 0.0f };
	}

	return matrix;
}

Matrix Matrix::CreateIdentityMatrix(size_t length)
{
	Matrix matrix = { };

	matrix.rowNum = length;
	matrix.columnNum = length;
	matrix.elements = new float*[length];

	for (int i = 0; i < length; ++i)
	{
		matrix.elements[i] = new float[length] { 0.0f };
		matrix.elements[i][i] = 1.0f;
	}

	return matrix;
}

size_t Matrix::GetRowNum() const
{
	return rowNum;
}

size_t Matrix::GetColumnNum() const
{
	return columnNum;
}

void Matrix::SetEpsilon(float epsilon)
{
	this->epsilon = epsilon;
}

float Matrix::GetEpsilon() const
{
	return epsilon;
}

void Matrix::GetDiagonalElements(float diagonalElements[]) const
{
	int min = rowNum < columnNum ? rowNum : columnNum;

	for (int i = 0; i < min; ++i)
	{
		diagonalElements[i] = elements[i][i];
	}
}

Matrix Matrix::GetTransposeMatrix() const
{
	Matrix matrix = CreateZeroMatrix(columnNum, rowNum);

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			matrix[j][i] = elements[i][j];
		}
	}

	return matrix;
}

Matrix Matrix::GetSubmatrix(size_t row, size_t column) const
{
	assert(row < rowNum && column < columnNum, "Cannot create a submatrix.");

	Matrix matrix = CreateZeroMatrix(rowNum - 1, columnNum - 1);

	int i = 0;
	int j = 0;
	int k = 0;
	int l = 0;

	for (i = 0; i < rowNum; ++i)
	{
		if (i == row)
		{
			continue;
		}

		for (j = 0, l = 0; j < columnNum; ++j)
		{
			if (j != column)
			{
				matrix[k][l] = elements[i][j];
				++l;
			}
		}

		++k;
	}

	return matrix;
}

float Matrix::GetDeterminant() const
{
	assert(IsSquareMatrix(), "Determinant does not exist");

	float total = 0.0f;
	float confactor;

	if (elements == nullptr)
	{
		return 0.0f;
	}
	else if (rowNum <= 1)
	{
		return **elements;
	}
	else
	{
		for (int i = 0; i < columnNum; ++i)
		{
			confactor = elements[0][i] * GetSubmatrix(0, i).GetDeterminant();
			confactor *= (i % 2 == 0) ? 1.0f : -1.0f;
			total += confactor;
		}

		return total;
	}
}

Matrix Matrix::GetAdjointMatrix() const
{
	assert(IsSquareMatrix(), "Determinant does not exist");

	Matrix adjointMatrix = CreateZeroMatrix(rowNum, columnNum);

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			adjointMatrix[i][j] = GetSubmatrix(i, j).GetDeterminant();
			adjointMatrix[i][j] *= ((i + j) % 2 == 0) ? 1.0f : -1.0f;
		}
	}

	return adjointMatrix.GetTransposeMatrix();
}

Matrix Matrix::GetInverseMatrix() const
{
	return GetAdjointMatrix() / GetDeterminant();
}

bool Matrix::IsZeroMatrix() const
{
	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			if (fabs(elements[i][j]) > epsilon)
			{
				return false;
			}
		}
	}

	return true;
}

bool Matrix::IsSymmetricMatrix() const
{
	if (!IsSquareMatrix())
	{
		return false;
	}

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = i + 1; j < columnNum; ++j)
		{
			if (fabs(elements[i][j] - elements[j][i]) > epsilon)
			{
				return false;
			}
		}
	}

	return true;
}

bool Matrix::IsSkewSymmetricMatrix() const
{
	if (!IsSquareMatrix())
	{
		return false;
	}

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = i + 1; j < columnNum; ++j)
		{
			if (fabs(elements[i][j] + elements[j][i]) > epsilon)
			{
				return false;
			}
		}
	}

	return true;
}

bool Matrix::IsSquareMatrix() const
{
	return rowNum == columnNum;
}

bool Matrix::IsDiagonalMatrix() const
{
	if (!IsSquareMatrix())
	{
		return false;
	}

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			if (fabs(elements[i][j]) > epsilon && i != j)
			{
				return false;
			}
		}
	}

	return true;
}

bool Matrix::IsIdentityMatrix() const
{
	if (!IsDiagonalMatrix())
	{
		return false;
	}

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			if (fabs(elements[i][j] - 1.0f) > epsilon)
			{
				return false;
			}
		}
	}

	return true;
}

Matrix& Matrix::operator=(const Matrix& matrix)
{
	assert(rowNum == matrix.rowNum && columnNum == matrix.columnNum, "Cannot perform operation due to different sizes of matrices.");

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			elements[i][j] = matrix[i][j];
		}
	}

	return *this;
}

Matrix Matrix::operator+(const Matrix& matrix) const
{
	assert(rowNum == matrix.rowNum && columnNum == matrix.columnNum, "Cannot perform operation due to different sizes of matrices.");

	Matrix ret = *this;
	
	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			ret.elements[i][j] += matrix.elements[i][j];
		}
	}

	return ret;
}

Matrix Matrix::operator-(const Matrix& matrix) const
{
	assert(rowNum == matrix.rowNum && columnNum == matrix.columnNum, "Cannot perform operation due to different sizes of matrices.");

	Matrix ret = *this;

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			ret.elements[i][j] -= matrix.elements[i][j];
		}
	}

	return ret;
}

Matrix Matrix::operator*(const Matrix& matrix) const
{
	assert(columnNum == matrix.rowNum, "Matrix multiplication does not satisfy the necessary conditions.");

	Matrix ret = CreateZeroMatrix(rowNum, matrix.columnNum);

	int i;
	int j;
	int k;
	int curRow = 0;
	int curColumn = 0;

	for (i = 0; i < rowNum; ++i)
	{
		curColumn = 0;

		for (j = 0; j < matrix.columnNum; ++j)
		{
			for (k = 0; k < columnNum; ++k)
			{
				ret[curRow][curColumn] += (elements[curRow][k] * matrix[k][curColumn]);
			}

			++curColumn;
		}

		++curRow;
	}

	return ret;
}

Matrix Matrix::operator/(const Matrix& matrix) const
{
	assert(columnNum == matrix.rowNum, "Matrix multiplication does not satisfy the necessary conditions.");

	return *this * matrix.GetInverseMatrix();
}

Matrix Matrix::operator*(float value) const
{
	Matrix ret = *this;

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			ret.elements[i][j] *= value;
		}
	}

	return ret;
}

Matrix Matrix::operator/(float value) const
{
	assert(value != 0.0f, "Zero division");

	Matrix ret = *this;

	int i;
	int j;

	for (i = 0; i < rowNum; ++i)
	{
		for (j = 0; j < columnNum; ++j)
		{
			ret.elements[i][j] /= value;
		}
	}

	return ret;
}

Matrix operator*(float value, const Matrix& matrix)
{
	Matrix ret = matrix;

	int i;
	int j;

	for (i = 0; i < ret.GetRowNum(); ++i)
	{
		for (j = 0; j < ret.GetColumnNum(); ++j)
		{
			ret.elements[i][j] *= value;
		}
	}

	return ret;
}

Matrix operator/(float value, const Matrix& matrix)
{
	assert(value != 0.0f, "Zero division");

	return matrix.GetInverseMatrix() * value;
}

Matrix& Matrix::operator+=(const Matrix& matrix)
{
	*this = *this + matrix;

	return *this;
}

Matrix& Matrix::operator-=(const Matrix& matrix)
{
	*this = *this - matrix;

	return *this;
}

Matrix& Matrix::operator*=(const Matrix& matrix)
{
	*this = *this * matrix;

	return *this;
}

Matrix& Matrix::operator/=(const Matrix& matrix)
{
	*this = *this / matrix;

	return *this;
}

Matrix& Matrix::operator*=(float value)
{
	*this = *this * value;

	return *this;
}

Matrix& Matrix::operator/=(float value)
{
	*this = *this / value;

	return *this;
}

float* Matrix::operator[](int i) const
{
	assert(i < rowNum, "Out of range.");

	return elements[i];
}