从头开始实现一个线性代数库：算法实现篇

最近重新学习了一下线性代数，读完了《程序猿的数学 3 线性代数》，为深刻理解内容，打算实现一个矩阵库，并写 Python 模块。计划涉及如下部分：

~~向量内积~~
~~矩阵乘法~~
~~LU 分解~~
~~线性方程组求解~~
~~行列式计算~~
~~逆矩阵~~
~~实对称矩阵 Jocobi 方法（求解特征值）~~
QR 分解

目前项目地址：https://github.com/netcan/LinAlg

关于 Python 模块的实现，可以参考我的这一篇文章：从头开始实现一个线性代数库：Python 模块篇

实现

内部统一采用 double 数据类型。

LU 分解

LU 分解将矩阵 A 分解为两个矩阵的乘积，其中为下三角矩阵，而为上三角矩阵，的形式如下，实现的时候采用一个矩阵来存储的结果。

LU 分解可用来方便的求解线性方程组、行列式的值以及逆矩阵。

分解方法：

将 L, U 改写成分块矩阵的形式：有

对，令

则即

对进行和一样的操作，最终得到零矩阵分解完毕可以得到, 。

核心代码如下：

Matrix Matrix::LUdecomp() const {
	Matrix LU(*this); // 最终结果存放到一个矩阵中
	size_t s = std::min(getRowSize(), getColSize());
	for(size_t k = 0; k < s; ++k) {
		double x = 1.0 / LU.data[k][k];
		// l1 = 1/a_{11} * ...
		for(size_t i = k + 1; i < getRowSize(); ++i)
			LU.data[i][k] *= x;
		// A -= lu^t,
		for(size_t i = k + 1; i < getRowSize(); ++i)
			for(size_t j = k + 1; j < getColSize(); ++j)
				LU.data[i][j] -= LU.data[i][k] * LU.data[k][j];
	}
	return LU;
}

若分解过程中出现除数为零，这需要添加一个置换矩阵，分解得到的形式，详见Matrix::PLUdecomp()。

行列式计算

对进行分解后，有

值即为对角线上元素的乘积。核心代码如下：

double Matrix::det() const {
	assert(getRowSize() == getColSize());
	auto PLU = PLUdecomp();
	const auto &P = PLU.first;
	const auto &LU = PLU.second;

	bool sign = false; // 负号
	for(size_t i = 0; i < PLU.first.size(); ++i)
		for(size_t j = i + 1; j < PLU.first.size(); ++j)
			if(PLU.first[i] > PLU.first[j]) sign = !sign;
	double ret = sign?-1.0:1.0;
	// det(A) = det(PLU) = det(P)det(L)det(U) = det(P)det(U)
	for(size_t i = 0; i < LU.data.size(); ++i)
		ret *= LU.data[P[i]][i];
	return ret;
}

求解线性方程组

对进行分解后，有

先用解出，在用解出。在解时，从头开始计算；而在解时，从尾开始计算。核心代码如下：

Vector Matrix::LUsolve_L(const Matrix &LU, const Vector &y) const {
	// 求出 Lz = y 的 z
	Vector z(y.getSize());
	for(size_t i = 0; i < z.getSize(); ++i) {
		z.data[i] = y.data[i];
		for (size_t j = 0; j < i; ++j)
			z.data[i] -= LU.data[i][j] * z.data[j];
	}
	return z;
}

Vector Matrix::LUsolve_U(const Matrix &LU, const Vector &z) const {
	// 求出 Ux = z 的 x
	Vector x(z.getSize());
	for(int i = z.getSize() - 1; i >=0; --i) {
		x.data[i] = z.data[i];
		for(int j = i + 1; j < z.getSize(); ++j)
			x.data[i] -= LU.data[i][j] * x.data[j];
		x.data[i] /= LU.data[i][i];
	}
	return x;
}

Vector Matrix::LUsolve(const Vector &y) const {
	// Ax = y
	assert(y.getSize() == getRowSize());
	assert(! isZero(det()) && "det is equal 0");
	Matrix LU = std::move(LUdecomp());
	Vector x(y.getSize());
	x = std::move(LUsolve_L(LU, y));
	x = std::move(LUsolve_U(LU, x));
	return x;
}

求解逆矩阵

对于，有 $展开得$

利用前面求解线性方程组的方法，得到。核心代码如下：

Matrix Matrix::inv() const {
	// A(x_1, ..., x_n) = (e_1, ..., e_n)
	assert(! isZero(det()) && "matrix is not invertible");
	size_t n = getRowSize();
	Matrix LU = std::move(LUdecomp());
	Matrix ret(n, n);
	Vector e(n), x(n);
	e[0] = 1.0;
	for(size_t j = 0; j < n; ++j) {
		if(j >= 1) std::swap(e[j], e[j-1]);
		// 解出 Ax = e
		x = std::move(LUsolve_L(LU, e));
		x = std::move(LUsolve_U(LU, x));
		for(size_t i = 0; i < n; ++i)
			ret.data[i][j] = x[i];
	}
	return ret;
}

利用 Jacobi 算法求解实对称矩阵特征值

对于 实对称 矩阵，将其对角化，有，这里的为对角矩阵，对角线上的值即为特征值，为特征向量组成的矩阵。

而解实对角矩阵的特征值算法有 Jacobi，主要通过不断平面旋转变化，使得矩阵朝着的方向变换。

定义如下平面旋转矩阵（正交矩阵）： $第列第列第行第行$

易知，所以有。通过正交矩阵将对称矩阵不断旋转变换进行对角化，。对左乘后，两行将发生变化，而右乘后，两列发生变化。

将 p, q 两行两列的将四个交叉点的结果写出来，由于只有这四个交叉点同时受到行列影响，所以它们的公式略有不同（由于是对称矩阵有）：

经过反复旋转变换后，矩阵依旧保持着对称性，依次选择，令旋转变化为得到，在这个作用下更新其他的行、列值。简而言之，每次将两行两列相交的四个交叉点中的斜对角交叉点变为 0，最终趋向于对角矩阵。

通过得到的值，有：

设，当不在增大时，可视为近似于对角矩阵。核心代码如下：

void Matrix::R(double cosphi, size_t p, size_t q, bool T) {
	// 通过平面旋转进行相似变换
	double sinphi = sqrt(1 - cosphi * cosphi);
	size_t n = getRowSize();
	for(size_t k = 0; k < n; ++k) {
		double  ap =  data[T?p:k][T?k:p] * cosphi + data[T?q:k][T?k:q] * sinphi,
				aq = -data[T?p:k][T?k:p] * sinphi + data[T?q:k][T?k:q] * cosphi;
		data[T?p:k][T?k:p] = ap; data[T?q:k][T?k:q] = aq;
	}
	return;
}
vector<double> Matrix::Jacobi() const {
	assert(getRowSize() == getColSize());
	size_t n = getRowSize();
	// 判断是否对称
	for(size_t i = 0; i < n; ++i)
		for(size_t j = 0; j < n; ++j)
			assert(Eq(data[i][j], data[j][i]));
	vector<double> eigenvals(n);
	Matrix D(*this); // 对角化
	// g(A) = \sum a[i][i]^2
	auto g = [](const Matrix &A)->double {
		double ret = 0.0;
		for(size_t i = 0; i < A.getRowSize(); ++i)
			ret += A.data[i][i] * A.data[i][i];
		return ret;
	};
	// f(A) -> 0 表明近似于对角化矩阵
	auto f = [](const Matrix &A)->double {
		double ret = 0.0;
		for(size_t i = 0; i < A.getRowSize(); ++i)
			for(size_t j = 0; j < A.getRowSize(); ++j)
				if(i!=j) ret += A.data[i][j] * A.data[i][j];
		return ret;
	};

	double cur_g = g(D), last_g = 0.0;

	while (! Eq(cur_g, last_g)) { // 当 g(D)不变时，对角化完成
		for (size_t p = 0; p < n; ++p) { // 逐个处理 a_{pq}
			for (size_t q = p + 1; q < n; ++q) {
				if (isZero(D.data[p][q])) continue;
				double  tan2phi = 2 * D.data[p][q] / (D.data[p][p] - D.data[q][q]),
						cosphi = sqrt(0.5 * (1.0 + 1.0/sqrt(1+tan2phi * tan2phi)));
				// D = R^T A R，旋转变换
				D.R(cosphi, p, q, true);
				D.R(cosphi, p, q, false);
			}
		}
		last_g = cur_g;
		cur_g = g(D);
	}
	for(size_t i = 0; i < n; ++i) eigenvals[i] = D.data[i][i];
	return eigenvals;
}