ROOT
ROOT
文章目录
  1. 实现
  2. LU 分解
    1. 分解方法:
  3. 行列式计算
  4. 求解线性方程组
  5. 求解逆矩阵
  6. 利用 Jacobi 算法求解实对称矩阵特征值

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

最近重新学习了一下线性代数,读完了《程序猿的数学 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;
}
支持一下
扫一扫,支持Netcan
  • 微信扫一扫
  • 支付宝扫一扫