- 1关系型数据库架构介绍及主流应用场景_无共享master和有共享master区别
- 2oauth2+JWT实现oauth2服务
- 3Docker- chapter 1
- 4【deeplearing.ai】 Week 16 RNN进阶以及Attention模型_imagecapture模型
- 5LATEX——加入图片_texstudio插入图片
- 6实习工作经历：代码在本地明明可以跑通，怎么放到服务器上就不行了呢？_我写的jq代码为什么放在服务器不能直接运行呢?
- 7零基础学编程从哪里入手，编程实例分享，配件进出库管理系统软件
- 8判断素数/质数的快速算法_prime/number/list
- 9遥感图像去雾文章解读_haze detection and removal in remotely sensed mult
- 10表白专用动画_var b = document.body; var c = document.getelement
六种语言求解特征值和特征向量（2022.4.29）_求矩阵特征值的程序代码
赞
踩
特征值和特征向量求解（2022.4.29）
引言
在使用层次分析法（AHP）来确定指标权重时，需要计算二维矩阵的特征值和特征向量，同时特征值和特征向量在计算机图像处理、物理学等领域都要十分重要的应用，下面首先简单介绍其定义和性质，之后利用六种编程语言进行实例求解。
1、矩阵的特征值和特征向量
A
=
[
a
11
a
12
.
.
.
a
1
n
a
21
a
22
.
.
.
a
2
n
.
.
.
.
.
.
.
.
.
.
.
.
a
n
1
a
n
2
.
.
.
a
n
n
]
A=\left [ a11a12...a1na21a22...a2n............an1an2...ann \right ]
A=⎣⎢⎢⎡a11a21...an1a12a22...an2............a1na2n...ann⎦⎥⎥⎤ 对于n阶方阵A，若存在数λ和一个n维的非零列向量
x
⃗
\vec{x}
x
，满足
A
x
⃗
=
λ
x
⃗
A\vec{x}=\lambda \vec{x}
Ax
=λx
，则称
x
⃗
\vec{x}
x
为矩阵A的对应λ特征值的特征向量，E为单位矩阵（主对角线元素均为1，其余为0）。
一方面，
(
λ
E
−
A
)
x
⃗
=
0
(\lambda E -A) \vec{x}=0
(λE−A)x
=0是包含n个未知数n个方程的齐次线性方程组，相当于存在不全为0的n个数使得矩阵
(
λ
E
−
A
)
(\lambda E -A)
(λE−A)的列向量线性相关，那么矩阵
(
λ
E
−
A
)
(\lambda E -A)
(λE−A)非满秩，即
r
a
n
k
(
λ
E
−
A
)
<
n
rank(\lambda E -A)<n
rank(λE−A)<n 另一方面，线性代数中齐次线性方程组存在非零解的充要条件是：系数矩阵行列式的值为零。因此得到下面的系数行列式：
∣
λ
E
−
A
∣
=
0
|\lambda E -A|=0
∣λE−A∣=0
∣
λ
E
−
A
∣
=
[
λ
−
a
11
−
a
12
.
.
.
−
a
1
n
−
a
21
λ
−
a
22
.
.
.
−
a
2
n
.
.
.
.
.
.
.
.
.
.
.
.
−
a
n
1
−
a
n
2
.
.
.
λ
−
a
n
n
]
=
λ
n
+
α
1
λ
n
−
1
+
α
2
λ
n
−
2
+
.
.
.
+
α
n
−
1
λ
+
α
n
=
0
|\lambda E -A|=\left [ λ−a11−a12...−a1n−a21λ−a22...−a2n............−an1−an2...λ−ann \right ]=\lambda^{n}+\alpha_{1}\lambda^{n-1}+\alpha_{2}\lambda^{n-2}+...+\alpha_{n-1}\lambda+\alpha_{n}=0
∣λE−A∣=⎣⎢⎢⎡λ−a11−a21...−an1−a12λ−a22...−an2............−a1n−a2n...λ−ann⎦⎥⎥⎤=λn+α1λn−1+α2λn−2+...+αn−1λ+αn=0 上式也称为方阵A的特征多项式，此多项式是一个关于参数λ的n次多项式。在复数域内（实数和虚数），n次代数方程有且仅有n个解。那么特征多项式的n个解可设为
λ
1
、
λ
2
、
⋯
、
λ
n
\lambda_{1}、\lambda_{2}、\cdots、\lambda_{n}
λ1、λ2、⋯、λn，这些解也称特征根，它们具有如下性质：
(
λ
−
λ
1
)
(
λ
−
λ
2
)
⋯
(
λ
−
λ
n
)
=
0
(\lambda -\lambda_{1})(\lambda -\lambda_{2})\cdots (\lambda -\lambda_{n})=0
(λ−λ1)(λ−λ2)⋯(λ−λn)=0
λ
1
+
λ
2
+
.
.
.
+
λ
n
=
∑
i
=
1
n
a
i
i
\lambda_{1}+\lambda_{2}+...+\lambda_{n}=\sum_{i=1}^{n}a_{ii}
λ1+λ2+...+λn=i=1∑naii
λ
1
λ
2
⋯
λ
n
=
∣
A
∣
\lambda_{1}\lambda_{2}\cdots\lambda_{n}=|A|
λ1λ2⋯λn=∣A∣
2、六种编程语言求解特征值和特征向量代码
为了求解如下方阵A的特征值和特征向量，下面分别利用Matlab、Python、C++、Java、C#、和R这六种编程语言进行代码实现，同时对比测试结果。
A
=
[
3
4
5
10
3
7
1
7
9
2
6
2
4
8
5
2
4
7
1
9
4
12
78
32
65
31
23
16
95
74
36
21
28
49
58
12
18
16
17
34
52
29
2
7
99
31
57
90
83
]
A=\left [ 3451037179262485247194127832653123169574362128495812181617345229279931579083 \right ]
A=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡3751295122492787418752432361699106765211731321312834577492349529018416582983⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
2.1 Matlab代码及运行结果
作为工业界和学术界的神器，MathWorks旗下的Matlab备受科研人员、工程师和开发者的喜爱，功能强大且有着“矩阵实验室”称号，在数学运算、图像处理、信号仿真方面有着其他软件无法比拟的计算优势，尽管有时其效率不是太高，但其拥有着毋庸置疑的准确率和精度。
n = 7; A=[3,4,5,10,3,7,1; 7,9,2,6,2,4,8; 5,2,4,7,1,9,4; 12,78,32,65,31,23,16; 95,74,36,21,28,49,58; 12,18,16,17,34,52,29; 2,7,99,31,57,90,83]; [eigenvector,x] = eig(A); lambda = diag(x); fprintf('特征值：'); for i=1:n if i == n fprintf(['λ',num2str(i),'= ',num2str(lambda(i)),'\n']); else fprintf([' ','λ',num2str(i),'= ',num2str(lambda(i)),' ']); end end for i=1:n if i == n fprintf(['特征向量',num2str(i),'\n']); else fprintf([' ','特征向量',num2str(i),' ']); end end disp(eigenvector);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
Matlab代码运行结果如下图所示：
2.2 Python代码及运行结果（Numpy）
这里利用Python 3.6和numpy依赖库来对矩阵进行特征值和特征向量运算，具体代码如下。
import numpy as np n = 7 A = np.array([[3,4,5,10,3,7,1], [7,9,2,6,2,4,8], [5,2,4,7,1,9,4], [12,78,32,65,31,23,16], [95,74,36,21,28,49,58], [12,18,16,17,34,52,29], [2,7,99,31,57,90,83]]) print('方阵A：\n{}'.format(A)) eigenvalue, eigenvector = np.linalg.eig(A) print('特征值', end='') for i in range(0, n): if i == n-1: print('λ{}={:.4f}'.format(i + 1, eigenvalue[i])) else: print('λ{}={:.4f}'.format(i+1, eigenvalue[i]), end=', ') for i in range(0, n): if i == n-1: print(' 特征向量{} '.format(i + 1)) else: print(' 特征向量{} '.format(i + 1), end=' ') for i in range(0, n): for j in range(0, n): if j == n-1: print('{:.4f} '.format(eigenvector[i][j])) else: print('{:.4f} '.format(eigenvector[i][j]), end=' ')
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
Python代码运行结果如下图所示：
2.3 C++代码及运行结果(Eigen)
Eigen作为线性代数相关计算的C++模板库广受开发人员的喜爱，下面使用Eigen库来计算矩阵特征值和特征向量，下载eigen-3.4.0后，在VS 2015中将其目录下的Eigen文件夹添加到项目属性的库目录中即可使用。
#include <iostream> #include <iomanip> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() { const int n = 7; Matrix<double, n, n> A = MatrixXd::Random(n, n); double a[] = { 3,4,5,10,3,7,1, 7,9,2,6,2,4,8, 5,2,4,7,1,9,4, 12,78,32,65,31,23,16, 95,74,36,21,28,49,58, 12,18,16,17,34,52,29, 2,7,99,31,57,90,83}; A = Map<Matrix<double, n, n, RowMajor> >(a); std::cout << "方阵A：" << std::endl; std::cout << A << std::endl; EigenSolver<Matrix<double, n, n>> eigensolver(A); MatrixXcd evecs = eigensolver.eigenvectors();//矩阵特征向量complex复数矩阵 MatrixXcd evals = eigensolver.eigenvalues(); //cout << evals << endl; //cout << evecs << endl; std::cout << "特征值："; for (int i = 0; i < n; i++) if(evals(i).imag() >= 0) cout << "λ" << i + 1 << "="<< setprecision(4) <<evals(i).real() <<"+"<< setprecision(4)<< evals(i).imag()<< "i "; else cout << "λ" << i + 1 << "="<< setprecision(4) << evals(i).real() << setprecision(4) << evals(i).imag() << "i "; std::cout << std::endl; for(int i=0;i<n;i++) std::cout << "特征向量" << (i+1) <<" "; std::cout << std::endl; for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { if(evecs(i, j).imag() >= 0) std::cout << setprecision(4) << evecs(i, j).real() <<"+"<< setprecision(4) << evecs(i, j).imag()<<"i "; else std::cout << setprecision(4) << evecs(i, j).real() << setprecision(4) << evecs(i, j).imag() << "i "; } std::cout << std::endl; } return 0; }
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
C++代码运行结果如下图所示：
2.4 Java代码及运行结果(ujmp)
这里使用第三方jar包：ujmp-core-0.3.0.jar来实现Java代码，jar包下载地址为ujmp-core-0.3.0.jar，之后在Java项目中需要将其添加到Build Path中，才能在代码中导入。
package com.test; import java.text.DecimalFormat; import org.ujmp.core.DenseMatrix; import org.ujmp.core.Matrix; import org.ujmp.core.doublematrix.impl.ArrayDenseDoubleMatrix2D; public class EigenValue_EigenVector { public static void main(String[] args) { int n = 7; double[][] a = {{3,4,5,10,3,7,1}, {7,9,2,6,2,4,8}, {5,2,4,7,1,9,4}, {12,78,32,65,31,23,16}, {95,74,36,21,28,49,58}, {12,18,16,17,34,52,29}, {2,7,99,31,57,90,83}}; Matrix A = DenseMatrix.Factory.importFromArray(a); System.out.print("A:"); System.out.println(A); Matrix[] eigenValueDecomposition = A.eig(); Matrix eigen_vector = eigenValueDecomposition[0]; Matrix eigen_value = eigenValueDecomposition[1]; double[][] eig_val= eigen_value.toDoubleArray(); double[][] eig_vec= eigen_vector.toDoubleArray(); DecimalFormat decimalFormat = new DecimalFormat("###################.####"); System.out.print("特征值"); for(int i=0;i<n;i++) { String val = String.valueOf(decimalFormat.format(eig_val[i][i])); for(int j=0;j<n;j++) { if(j !=i && eig_val[j][i] != 0) if(eig_val[j][i] >= 0) val += "+"+String.valueOf(decimalFormat.format(eig_val[j][i]))+"i"; else val += String.valueOf(decimalFormat.format(eig_val[j][i]))+"i"; } System.out.print("λ"+Integer.toString(i+1)+"= "+val+" "); } System.out.println(); for(int i=0;i<n;i++) { System.out.print(" 特征向量"+Integer.toString(i+1)+" "); } System.out.println(); for(int i=0;i<n;i++) { for(int j=0;j<n;j++) { System.out.print("\t"+String.valueOf(decimalFormat.format(eig_vec[i][j]))+"\t"); } System.out.println(); } } }
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
Java代码运行结果如下图所示：
2.5 C#代码及运行结果(MathNet)
这里利用了适用于C#语言的经典数学dll库:MathNet.Numerics.dll，在编写代码时需要将其添加到引用当中，以便在代码中调用其接口。
using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Numerics; using MathNet.Numerics; using MathNet.Numerics.LinearAlgebra; using MathNet.Numerics.LinearAlgebra.Double;//必须 using MathNet.Numerics.LinearAlgebra.Complex; using MathNet.Numerics.LinearAlgebra.Factorization; namespace SolveMatrixEigen { class Program { public static void Main(string[] args) { Console.WriteLine("Hello，World!"); const int n = 7; double[,] a = new double[n, n]{{3,4,5,10,3,7,1}, {7,9,2,6,2,4,8}, {5,2,4,7,1,9,4}, {12,78,32,65,31,23,16}, {95,74,36,21,28,49,58}, {12,18,16,17,34,52,29}, {2,7,99,31,57,90,83}}; MathNet.Numerics.LinearAlgebra.Matrix<double> A = MathNet.Numerics.LinearAlgebra.Double.DenseMatrix.OfArray(a); MathNet.Numerics.LinearAlgebra.Factorization.Evd<double> eigen = A.Evd(); Vector<Complex> Eigen_Value = eigen.EigenValues;//矩阵的特征值 MathNet.Numerics.LinearAlgebra.Matrix<double> Eigen_Vector = eigen.EigenVectors;//矩阵的特征向量 Console.Write("特征值："); for (int i = 0; i < n; i++) { string eigenvalue_i = ""; if(Eigen_Value[i].Imaginary >= 0) eigenvalue_i = String.Format("{0:F4}+{1:F4}i", Eigen_Value[i].Real, Eigen_Value[i].Imaginary); else eigenvalue_i = String.Format("{0:F4}{1:F4}i", Eigen_Value[i].Real, Eigen_Value[i].Imaginary); Console.Write("λ" + (i + 1).ToString() + "=" + eigenvalue_i+" "); } Console.WriteLine(); for (int i = 0; i < n; i++) Console.Write(" 特征向量{0} ", i + 1); Console.WriteLine(); for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) Console.Write(String.Format(" {0:F4} ",Eigen_Vector[i,j])); Console.WriteLine(); } } } }
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
C#代码运行结果如下图所示：
2.6 R代码及运行结果
Windows下的R安装包下载，安装成功后出现RGui启动程序，进入后可直接编写脚本代码。
a <- c(3,4,5,10,3,7,1,7,9,2,6,2,4,8,5,2,4,7,1,9,4,12,78,32,65,31,23,16,95,74,36,21,28,49,58,12,18,16,17,34,52,29,2,7,99,31,57,90,83)
A <- matrix(a,7,7)
print("方阵A")
print(A)
ev <- eigen(A)
print("特征值：")
print(ev$values)
print("特征向量：")
print(ev$vectors)
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
R代码运行结果如下图所示：
3、总结
总的来说，Matlab、Python和R同为脚本语言，可解释性强，易于执行，易于验证，直观高效；而Java、C#和C++都为面向对象的高级编程语言，需要依托主函数才能执行代码，效率较高，计算速度快，但往往需要额外依托第三方库或者自己开发代码来执行矩阵运算。
同时实际测试发现：六种编程语言对特征值的计算都十分准确，但特征向量结果不同，具体不同之处在于：Matlab和Python计算的特征向量结果相同，为复数；Java和C#计算的特征向量结果相同，是实数；C++计算的特征向量结果为复数，但与Matlab和Python的计算结果略有差异；R处理结果相差较大。最终利用Matlab软件验证特征值和特征向量的结果如下图所示。因此，无论利用何种编程语言进行计算，如果对计算结果的精度要求较高，建议与专业软件进行结果比对和相互验证，从而保障计算的准确性，这也告诉我们第三方库也无法保证算法的高度准确，因为它们所采用的计算方式可能有所差异。
 |  |
 |  |
 |  |
 |
