稀疏张量 分解 python 稀疏求解算法

已知线性方程组的矩阵表示形式为：
$Ax=b$
当矩阵$A$为低阶稠密矩阵时，通常可以采用直接法-矩阵分解，将矩阵直接进行分解，然后求解；
当矩阵为高阶稀疏矩阵时（含有较多0元素），通常采用迭代法，如雅克比（Jacobi）迭代法、高斯-赛德尔（Gauss-Seidel）迭代法、超松弛（SOR）迭代法等，对于迭代法，将会面临收敛性问题。即随着迭代次数的增加，误差并不会较小，反而不断增加。

本文分别对雅克比迭代法、高斯-赛德尔迭代法和超松弛迭代法，用矩阵形式进行表示，并分析它们的收敛性。

全文定义

• 线性方程组：$Ax=b$
• 迭代公式：$x^{k+1}=Bx^{k}+g$
• 对于迭代公式$x^{k+1}=Bx^{k}+g$收敛的充分必要条件是跌倒矩阵$B$的谱半径$\rho(G)<1$
• 如果迭代矩阵的某种范数$||B||<1$，则迭代公式$x^{k+1}=Bx^{k}+g$收敛
• 对角占优的线性方程组$Ax=b$的雅克比迭代公式和高斯-赛德尔迭代公式都收敛。
• 若系数矩阵$A$是对称正定矩阵，则高斯-赛德尔迭代公式收敛

补充：

• 谱半径：矩阵$A$的谱半径$\rho(A)$就是矩阵$A$的特征值的模的最大值。设$A$是$n$阶方阵，$\lambda _i$是其特征值，则$\rho(A) = max\{|\lambda_i|\}$。
• 矩阵范数：

定义一个矩阵$A = \begin{bmatrix} -1&2&-3\\ 4&-6&6 \end{bmatrix}$

• 1范数：矩阵的每一列上的元素绝对值先求和，再从中取个最大的（列和最大），上述矩阵A的1范数先得到[5,8,9]，再取最大的最终结果就是：9。
• 2范数：矩阵$A^TA$的最大特征值开平方根，上述矩阵A的2范数得到的最大结果是：10.0623。
• 无穷范数：矩阵的每一行上的元素绝对值先求和，再从中取个最大的，（行和最大），上述矩阵A的1范数先得到[6;16],再取最大的最终结果就是：16。
• 核范数：矩阵的奇异值（将矩阵$SVD$分解）之和，这个范数可以用来低秩表示（因为最小化核范数，相当于最小化矩阵的秩–低秩），上述矩阵A的最终结果就是10.9287。
• $L_0$范数：矩阵的非0元素的个数，通常用它来表示稀疏，L0范数越小0元素越多，也就越稀疏，上述矩阵A最终结果就是：6。
• $L_1$范数：矩阵中的每个元素绝对值之和，它是$L_0$范数的最优凸近似，因此它也可以表示稀疏，上述矩阵A最终结果就是：22。
• $F$范数：矩阵的各个元素平方之和再开平方根，它通常也叫做矩阵的$L_2$范数，它的优点在它是一个凸函数，可以求导求解，易于计算，上述矩阵A最终结果就是10.0995。
• $L_{21}$范数：矩阵先以每一列为单位，求每一列的$F$范数（也可认为是向量的2范数），然后再将得到的结果求$L_1$范数（也可以认为是向量的1范数），很容易看出它是介于$L_1$和$L_2$之间的一种范数，上述矩阵A最终结果就是：17.1559。

1 雅克比迭代法

1.1 矩阵表示

雅克比迭代法，假设系数矩阵$A$非奇异，并且主对角线上的元素不等于0.
因此可将$A$矩阵进行分解：
$A = L+D+U$
其中
$L$为下三角矩阵，下三角各元素对应于系数矩阵$A$下三角上的元素，并且主对角线上的元素为0；
$D$为对角矩阵，主对角线上的元素即为系数矩阵$A$对角线上的元素。
$U$为上三角矩阵，上三角各元素对应于系数矩阵$A$上三角上的元素，并且主对角线上的元素为0；

通过推导，可得雅克比迭代法的迭代公式为：
$x^{k+1}=Bx^{k}+g\\B=(I-D^{-1}A)\\g=D^{-1}b\\x^{k+1} = (I-D^{-1}A)x^{k}+D^{-1}b$
其中：矩阵$B$称为雅克比迭代矩阵。

1.2 雅克比迭代法分量表示

雅克比迭代法的矩阵表示法，主要用来讨论起收敛性；实际计算中，要用雅克比迭代法公式的分量形式。

2 高斯-赛德尔迭代法

2.1 矩阵表示

$Ax=b=(L+D+u)x\\ Dx^{k+1}=-Lx^{k+1}-Ux^{k}+b$
因此高斯-赛德尔迭代法的矩阵表示为：
$x^{k+1}=Bx^{k}+g\\x^{k+1} = -(D+L)^{-1}Ux^{k}+(D+L)^{-1}b\\ B=-(D+L)^{-1}U\\ g=(D+L)^{-1}b$

2.2 高斯-赛德尔迭代法的分量表示：

它充分了利用了本次迭代计算的结果，即计算$x^{k+1}$的第$i$个分量时，由于它需要知道，此时其它分量的值，对于排在它之前的分量$x_j^{k+1}(j<i)$，由于在本次迭代过程中，已经计算出来本次迭代的结果，因此该分量$x_i^{k+1}$将会利用到$x_j^{k+1}(j<i)$；而对于排在它之后的分量$j>i$，因为本次迭代还未轮到它的计算，因此对于$x_i^{k+1}$的计算，只用利用第$j$分分量上一次迭代的结果$x_j^{k}$

3 逐次超松弛迭代法

$SOR$迭代法是高斯-赛德尔迭代法的一种修正，它通过$x^{k}$计算$x^{k+1}$的步骤和迭代公式为：

1. 利用高斯迭代法求得辅助量
2. 与$x_i^{k}$做加权平均

整合为一个公式即为：

其中：$\omega$为松弛因子。

3.1 $SOR$迭代法的矩阵表示

$x^{k+1} = Bx^{k}+g\\ x^{k+1} = (D-\omega L)^{-1}[(1-\omega)D+\omega U]x^{k}+\omega (D-\omega L)^{-1} b\\ B = (D-\omega L)^{-1}[(1-\omega)D+\omega U]\\g=\omega (D-\omega L)^{-1}$

• 当$\omega =1$时，即为高斯-赛德尔迭代法
• $SOR$迭代法每次迭代法的主要运算量为矩阵和向量的乘法
• 当$\omega>1$时，称为超松弛迭代法；当$\omega <1$时，称为低松弛迭代法。对于特定线性方程组的求解，合理选择松弛因子$\omega$的大小，将会使得$SOR$收敛大大加速。

$SOR$ 迭代法的定理：

• 若矩阵$A$为对称正定矩阵，当$0<\omega<2$时，$SOR$迭代法是收敛的
• 若矩阵$A$是严格的对角占优矩阵，当$0<\omega<1$时，，$SOR$迭代法是收敛的

4 迭代终止条件

或者

对于$x$的每一个分量$x_i$，它们的迭代次数$k$可能是不同的。因此不能通过矩阵形式直接求解。

5 代码实现

5.1 实现雅克比迭代法

#include<iostream>
#include<Eigen/Core>
#include<Eigen/Dense>
#include<chrono>

using namespace std;

typedef Eigen::Matrix<float, 6, 6>Matrix6f;
typedef Eigen::Matrix<float, 6, 1>Vector6f;
const double e = 0.000001; //最大误差
int main(int argc, char** argv)
{
	Matrix6f A;
	A << 4, -1, 0, -1, 0, 0,
		-1, 4, -1, 0, -1, 0,
		0, -1, 4, -1, 0, -1,
		-1, 0, -1, 4, -1, 0,
		0, -1, 0, -1, 4, -1,
		0, 0, -1, 0, -1, 4;
	Vector6f b;
	b << 0, 5, -2, 5, -2, 6;

	//cout << A << endl;
	//cout << b << endl;
	double er=1;
	//构造雅克比迭代公式
	Vector6f xk = Vector6f::Zero(); //将迭代初始值设为0
	int k = 0; //计算迭代的次数
	chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
	while (er>e)
	{
		er = 0;
		Vector6f  xk_1;
		Vector6f index;
		for (int i = 0; i < A.cols(); i++)
		{
			index = Vector6f::Zero();
			for (int j = 0; j < A.cols(); j++)
			{
				if (j != i)
				{
					index(i,0) = index(i,0) + A(i,j) * xk(j,0);
				}
			}
			xk_1(i,0) = (b(i,0) - index(i,0)) / A(i,i);
		}
		//计算误差
		er = (xk - xk_1).dot(xk - xk_1) / 6;
		xk = xk_1;
		k++;
	}
	chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
	//计算迭代法的时间
	chrono::duration<float> time1 = chrono::duration_cast<chrono::duration<float>>(t2 - t1);
	cout << "迭代的次数：";
	cout << k << endl;
	cout << "迭代的时间为：" << time1.count()<<"秒"<<endl;
	cout << "结果：" << endl;
	cout << xk;
	cout << endl;
	cout << "=====================";
	cout << endl;

	//与直接三角分级法相比
	Eigen::FullPivLU<Matrix6f>lu(A);

	chrono::steady_clock::time_point t3 = chrono::steady_clock::now();
	Vector6f lx = lu.solve(b);
	chrono::steady_clock::time_point t4 = chrono::steady_clock::now();
	chrono::duration<float> time2 = chrono::duration_cast<chrono::duration<float>>(t4 - t3);
	cout << "通过三角分解法，LU的计算时间为：" << time2.count();
	cout << endl;
	cout << "计算结果为：" << endl;
	cout << lx << endl;
	cout << endl;
	system("pause");
	return 0;
}

计算结果为：

迭代的次数：20
迭代的时间为：0.0014122秒
结果：
0.999615
 1.99895
0.999474
 1.99895
0.999474
 1.99923
=====================
通过三角分解法，LU的计算时间为：0.0001018
计算结果为：
1
2
1
2
1
2

5.2 实现高斯-赛德尔迭代法

//实现Gauss-Seidel迭代法
#include<iostream>
#include<Eigen/Core>
#include<Eigen/Dense>
#include<chrono>

using namespace std;

const double er = 0.0001;
typedef Eigen::Matrix<float, 6, 6> Matrix6f;
typedef Eigen::Matrix<float, 6, 1>Vector6f;
//Vector6f GS(const Matrix6f&, const Vector6f&); //高斯-赛德尔迭代函数声明

int main(int argc, char** argv)
{
	//构造数据
	Matrix6f A;
	Vector6f b;
	A << 4, -1, 0, -1, 0, 0,
		-1, 4, -1, 0, -1, 0,
		0, -1, 4, -1, 0, -1,
		-1, 0, -1, 4, -1, 0,
		0, -1, 0, -1, 4, -1,
		0, 0, -1, 0, -1, 4;
	b << 0, 5, -2, 5, -2, 6;
	
	//使用高斯-赛德尔迭代法计算线性方程组的解
	//定义初始解
	Vector6f x0 = Vector6f::Zero();
	Vector6f index = Vector6f::Zero();
	double e = 1;
	//计数迭代次数
	int k = 0;
	chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
	while (e > er)
	{
		e = 0;
		for (int i =0;i< A.cols(); i++)
		{
			index = Vector6f::Zero();
			for (int j = 0; j < A.cols(); j++)
			{
				if (j != i)
				{
					index(i, 0) += A(i, j)*x0(j, 0);
				}
			}
			x0(i, 0) = (b(i, 0) - index(i, 0)) / A(i, i);
		}
		//计算误差
		Vector6f new_b = A * x0;
		e = (new_b - b).dot(new_b - b) / 6;
		k += 1;
	}
	chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
	cout << "Gauss-Seidel迭代法和Jacobi迭代法" << endl;
	cout << "===================================" << endl;
	cout << "高斯-赛德尔迭代法" << endl;
	cout << "计算的结果为：" << endl;
	cout << x0 << endl;
	//cout << result << endl;
	cout << "Ax=" << endl;
	//cout << A * result << endl;
	cout << A * x0 << endl;
	cout << "迭代的次数为：" << k << endl;
	double time1 = chrono::duration_cast<chrono::duration<double>>(t2 - t1).count();
	cout << "迭代的时间为：" << time1<<endl;
	cout << "====================================" << endl;

	//Jacobi迭代法
	int kj = 0;
	Vector6f xj = Vector6f::Zero();
	Vector6f xj_1;
	Vector6f index1 = Vector6f::Zero();
	double ej = 1;
	chrono::steady_clock::time_point t3 = chrono::steady_clock::now();
	while (ej > er)
	{
		ej = 0;
		for (int i = 0; i < A.cols(); i++)
		{
			index1 = Vector6f::Zero();
			for (int j = 0; j < A.cols(); j++)
			{
				if (i != j)
				{
					index1(i, 0) += A(i, j)*xj(j, 0);
				}
			}
			//计算x_(k+1)
			xj_1(i, 0) = (b(i,0)-index1(i,0)) / A(i, i);
		}
		//计算两个解的差值
		ej = (xj - xj_1).dot(xj - xj_1) / 6;
		//每次迭代完成，将所有xj_1的解一次性复制给xj，用于下次求解
		xj = xj_1;
		xj_1 = Vector6f::Zero();
		kj += 1;
	}
	chrono::steady_clock::time_point t4 = chrono::steady_clock::now();
	double time2 = chrono::duration_cast<chrono::duration<double>>(t4 - t3).count();
	cout << "Jacobi方法：" << endl;
	cout << "计算的结果：" << endl;
	cout << xj << endl;
	cout << "Ax:" << endl;
	cout << A * xj << endl;
	cout << "迭代的次数：" << kj << endl;
	cout << "迭代的时间" << time2<<endl;
	cout << "==============================" << endl;
	system("pause");
	return 0;
}

Gauss-Seidel迭代法和Jacobi迭代法
===================================
高斯-赛德尔迭代法
计算的结果为：
0.993268
  1.9927
0.993837
 1.99436
 0.99536
  1.9973
Ax=
-0.0139883
   4.98833
  -2.00901
   4.99498
  -2.00292
         6
迭代的次数为：8
迭代的时间为：0.0005605
====================================
Jacobi方法：
计算的结果：
0.988895
 1.99242
0.984831
 1.99242
0.984831
 1.99445
Ax:
-0.0292492
    5.0111
  -2.03996
    5.0111
  -2.03996
   6.00813
迭代的次数：13
迭代的时间0.0007516
==============================
请按任意键继续. . .

5.3 实现SOR迭代法

#include<iostream>
#include<Eigen/Core>
#include<Eigen/Dense>
#include<chrono>

using namespace std;

typedef Eigen::Matrix<float, 6, 6> Matrix6f;
typedef Eigen::Matrix<float, 6, 1>Vector6f;
const double er = 0.0001;

int main(int argc, char** argv)
{
	//生成数据
	Matrix6f A;
	A << 4, -1, 0, -1, 0, 0,
		-1, 4, -1, 0, -1, 0,
		0, -1, 4, -1, 0, -1,
		-1, 0, -1, 4, -1, 0,
		0, -1, 0, -1, 4, -1,
		0, 0, -1, 0, -1, 4;
	Vector6f b;
	b << 0, 5, -2, 5, -2, 6;

	//定义误差
	double e = 1;
	//定义omega
	float omega;
	cout << "请输入w的值：";
	cin >> omega;
	
	//定义迭代初始值
	Vector6f xk = Vector6f::Zero();
	Vector6f xk_11;
	Vector6f xk_1 = Vector6f::Zero();
	//定义中间计算量
	Vector6f index = Vector6f::Zero();
	//定义迭代次数
	int k = 0;

	//SOR迭代
	chrono::steady_clock::time_point t1 = chrono::steady_clock::now();
	while (e > er)
	{
		e = 0;
		for (int i = 0; i < A.cols(); i++)
		{
			//保存迭代计算之前的值
			xk_11(i, 0) = xk(i, 0);
			index(i, 0) = 0;
			for (int j = 0; j < A.cols(); j++)
			{
				if (j != i)
				{
					index(i, 0) += A(i, j)*xk(j, 0);
				}
			}
			//计算x_{k+1}
			xk(i, 0) = (b(i, 0) - index(i, 0)) / A(i, i);
			//计算xk_1，每轮迭代完成的值都保存在vk_1向量中
			xk_1(i, 0) = (1 - omega)*xk_11(i, 0) + omega * xk(i, 0);
		}
		//更新每次迭代完成的信息
		//计算误差，前后的值分别保存在xk_11和xk_1中
		e = (xk_1 - xk_11).dot(xk_1 - xk_11) / 6;
		xk = xk_1;
		k += 1;
	}
	chrono::steady_clock::time_point t2 = chrono::steady_clock::now();
	double time = chrono::duration_cast<chrono::duration<double>>(t2 - t1).count();
	cout << "SOR迭代" << endl;
	cout << "解：" << endl;
	cout << xk;
	cout << "b" << endl;
	cout << A * xk << endl;
	cout << "迭代的次数：" << k << endl;
	cout << "迭代的时间：" << time << "秒" << endl;
	
	system("pause");
	return 0;
}

请输入w的值：1.2
SOR迭代
解：
0.996406
 1.99735
0.997802
 1.99783
0.998151
 1.99901b
-0.00956082
    4.99705
   -2.00298
    4.99897
   -2.00159
    6.00008
迭代的次数：7
迭代的时间：0.0004327秒
请按任意键继续. . .

请输入w的值：0.6
SOR迭代
解：
0.981769
 1.98007
0.983152
 1.98458
0.987272
  1.9926b
-0.0375741
   4.96807
  -2.02464
   4.98614
  -2.00816
   5.99996
迭代的次数：13
迭代的时间：0.0007727秒
请按任意键继续. . .